# Source code for naginterfaces.library.rnla

# -*- coding: utf-8 -*-
r"""
Module Summary
--------------
Interfaces for the NAG Mark 28.4 rnla Chapter.

rnla - Randomized Numerical Linear Algebra

This module covers linear algebra functions that make use of random projections to reduce problem dimension.
This area is referred to as RNLA (Randomized Numerical Linear Algebra).
The functions can be split into the following categories:

building blocks that are intended to be used as components in RNLA algorithms written by you;

RNLA algorithms for matrix factorization.

It is envisaged that users of the higher-level functions, such as matrix factorization, will have some background in linear algebra.
A common use case would be that you have tried solving your problem using a deterministic linear algebra function, e.g., an LAPACK function from submodule :mod:~naginterfaces.library.lapackeig, and are in need of a function that is more computationally efficient.

Users of the building block functions would be expected to have some familiarity with the RNLA literature, i.e., a higher level of expertise.

Functionality Index
-------------------

**Building blocks**

random projection with DCT (float): :meth:randproj_dct_real

**Matrix factorization**

SVD via row extraction (float): :meth:svd_rowext_real

For full information please refer to the NAG Library document

https://www.nag.com/numeric/nl/nagdoc_28.4/flhtml/f10/f10intro.html
"""

[docs]def svd_rowext_real(jobu, jobvt, a, k, statecomm, rtol_abs=-1.0, rtol_rel=-1.0):
r"""
svd_rowext_real computes the singular value decomposition (SVD) of a real :math:m\times n matrix :math:A, optionally computing the left and/or right singular vectors by using a randomised numerical linear algebra (RNLA) method.

.. _f10ca-py2-py-doc:

For full information please refer to the NAG Library document for f10ca

https://www.nag.com/numeric/nl/nagdoc_28.4/flhtml/f10/f10caf.html

.. _f10ca-py2-py-parameters:

**Parameters**
**jobu** : str, length 1
Specifies options for computing part of or none of the matrix :math:U.

:math:\mathrm{jobu} = \texttt{'S'}

The first :math:k columns of :math:U (the left singular vectors) are returned in the array :math:\mathrm{u}.

:math:\mathrm{jobu} = \texttt{'N'}

No columns of :math:U (no left singular vectors) are computed.

**jobvt** : str, length 1
Specifies options for computing part of or none of the matrix :math:V^\mathrm{T}.

:math:\mathrm{jobvt} = \texttt{'S'}

The first :math:k rows of :math:V^\mathrm{T} (the right singular vectors) are returned in the array :math:\mathrm{vt}.

:math:\mathrm{jobvt} = \texttt{'N'}

No rows of :math:V^\mathrm{T} (no right singular vectors) are computed.

**a** : float, array-like, shape :math:\left(m, n\right)
The :math:m\times n matrix :math:A.

**k** : int
:math:k, number of columns in random projection, :math:Y = A\Omega.

**statecomm** : dict, RNG communication object, modified in place
RNG communication structure.

This argument must have been initialized by a prior call to :meth:rand.init_repeat <naginterfaces.library.rand.init_repeat> or :meth:rand.init_nonrepeat <naginterfaces.library.rand.init_nonrepeat>.

**rtol_abs** : float, optional
The absolute tolerance, :math:\epsilon_a used in defining the threshold on estimating the rank of :math:A. If :math:\mathrm{rtol\_abs} < 0.0 then :math:0.0 is used unless :math:\mathrm{rtol\_rel}\leq 0.0 in which case :math:\sqrt{\text{machine precision}} is used.

**rtol_rel** : float, optional
The relative tolerance, :math:\epsilon_r used in defining the threshold on estimating the rank of :math:A. If :math:\mathrm{rtol\_rel} < 0.0 then :math:0.0 is used unless :math:\mathrm{rtol\_abs}\leq 0.0 in which case :math:\sqrt{\text{machine precision}} is used.

**Returns**
**s** : float, ndarray, shape :math:\left(\mathrm{r}\right)
The :math:\mathrm{r} largest singular values of :math:A in descending order.

**u** : None or float, ndarray, shape :math:\left(:, :\right)
If :math:\mathrm{jobu} = \texttt{'S'}, :math:\mathrm{u} contains the first :math:\mathrm{r} columns of :math:U (the left singular vectors, stored column-wise).

**vt** : None or float, ndarray, shape :math:\left(:, :\right)
If :math:\mathrm{jobvt} = \texttt{'S'}, :math:\mathrm{vt} contains the first :math:\mathrm{r} rows of :math:V^\mathrm{T} (the right singular vectors).

**r** : int
:math:r, contains estimated rank of array :math:\mathrm{a}.

.. _f10ca-py2-py-errors:

**Raises**
**NagValueError**
(errno :math:1)
On entry, :math:\mathrm{jobu} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{jobu} = \texttt{'S'} or :math:\texttt{'N'}.

(errno :math:2)
On entry, :math:\mathrm{jobvt} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{jobvt} = \texttt{'S'} or :math:\texttt{'N'}.

(errno :math:3)
On entry, :math:m = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:m > 0.

(errno :math:4)
On entry, :math:n = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:n > 0.

(errno :math:7)
On entry, :math:\mathrm{k} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:0 < \mathrm{k}\leq \mathrm{min}\left(m, n\right).

(errno :math:8)
On entry, :math:\mathrm{statecomm}\ ['state'] vector has been corrupted or not initialized.

(errno :math:9)
On entry, :math:\mathrm{jobu} = \langle\mathit{\boldsymbol{value}}\rangle and :math:\textit{ldu} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: if :math:\mathrm{jobu} = \texttt{'S'}, :math:\textit{ldu}\geq m.

(errno :math:10)
On entry, :math:\mathrm{jobvt} = \langle\mathit{\boldsymbol{value}}\rangle and :math:\textit{ldvt} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: if :math:\mathrm{jobvt} = \texttt{'S'}, :math:\textit{ldvt}\geq \mathrm{k}.

(errno :math:21)
On exit, :math:\mathrm{r} = \mathrm{k}, the rank of :math:A may be larger than :math:\mathrm{r}.

Increase :math:\mathrm{k} to obtain a more accurate rank estimate.

Smallest diagonal element of :math:R, from :math:QR of :math:Y, :math:=:math:\langle\mathit{\boldsymbol{value}}\rangle.

Tolerance used to determine rank :math:=:math:\langle\mathit{\boldsymbol{value}}\rangle.

(errno :math:22)
:math:A has effective rank of zero.

First diagonal element of :math:R, from :math:QR of :math:Y, :math:=:math:\langle\mathit{\boldsymbol{value}}\rangle.

Tolerance used to determine rank :math:=:math:\langle\mathit{\boldsymbol{value}}\rangle.

.. _f10ca-py2-py-notes:

**Notes**
The SVD is written as

.. math::
A = U\Sigma V^\mathrm{T}\text{,}

where :math:\Sigma is an :math:m\times n matrix which is zero except for its :math:\mathrm{min}\left(m, n\right) diagonal elements, :math:U is an :math:m\times m orthogonal matrix, and :math:V is an :math:n\times n orthogonal matrix.
The diagonal elements of :math:\Sigma are the singular values of :math:A; they are real and non-negative, and are returned in descending order.
The first :math:\mathrm{min}\left(m, n\right) columns of :math:U and :math:V are the left and right singular vectors of :math:A.

Note that the function returns :math:V^\mathrm{T}, not :math:V.

If the rank of :math:A is :math:r < \mathrm{min}\left(m, n\right), then :math:\sigma has :math:r nonzero elements, and only :math:r columns of :math:U and :math:V are well-defined.
In this case we can reduce :math:\sigma to an :math:r\times r matrix, :math:U to an :math:m\times r matrix and :math:V^\mathrm{T} to an :math:r\times n matrix.

svd_rowext_real is designed for efficiently computing the SVD in the case :math:r < \mathrm{min}\left(m, n\right).
The input argument :math:k should be greater than :math:r by a small oversampling parameter, :math:\delta, such that :math:k = r+\delta.
A reasonable value for :math:\delta, to compute the SVD to within machine precision, is :math:\delta = 5.
The value of :math:\delta should not vary based on :math:m or :math:n.
If :math:r is not known then the function can be used iteratively to refine the estimate and accuracy of the computed SVD; using a larger value of :math:\delta than necessary increases the computational cost of the function.

As a by-product of computing the SVD, the function estimates :math:r.

If the input argument :math:k is less than :math:r the accuracy depends on the :math:\left(k+1\right)\ th singular value, :math:\sigma_{{k+1}}.
See Accuracy <https://www.nag.com/numeric/nl/nagdoc_28.4/flhtml/f10/f10caf.html#accuracy>__ for more details.

A call to svd_rowext_real consists of the following:

(1) A random projection is applied, :math:Y = A\Omega, where :math:\Omega is an :math:n\times k matrix. (Note that the product :math:A\Omega is computed using a Fast Fourier Transform, so can be computed in :math:\mathrm{O}\left(mn\log\left(k\right)\right) time.) See :meth:randproj_dct_real for more details on the random projection.

(#) A pivoted :math:QR decomposition of :math:Y is calculated (see :meth:lapackeig.dgeqp3 <naginterfaces.library.lapackeig.dgeqp3> for more details). The rank estimate is then :math:r such that, on the diagonal of :math:R,

.. math::
\left\lvert R_{{rr}}\right\rvert > \epsilon_a+\epsilon_r\times R_{{11}}\text{,}

where :math:\epsilon_a and :math:\epsilon_r are the absolute and relative error tolerances, respectively, and :math:r is the largest diagonal index for which the above relation holds.

(#) Obtain the SVD from the :math:QR decomposition of :math:Y (or, depending on the rank, an approximation to the SVD) of :math:A. This is referred to as row extraction.

Further details of the randomized SVD procedure can be found in Sections 4 and 5 of Halko et al. (2011).

.. _f10ca-py2-py-references:

**References**
Anderson, E, Bai, Z, Bischof, C, Blackford, S, Demmel, J, Dongarra, J J, Du Croz, J J, Greenbaum, A, Hammarling, S, McKenney, A and Sorensen, D, 1999, LAPACK Users' Guide, (3rd Edition), SIAM, Philadelphia, https://www.netlib.org/lapack/lug

Golub, G H and Van Loan, C F, 1996, Matrix Computations, (3rd Edition), Johns Hopkins University Press, Baltimore

Halko, N, 2012, Randomized methods for computing low-rank approximations of matrices, PhD thesis

Halko, N, Martinsson, P G and Tropp, J A, 2011, Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions, SIAM Rev. (53(2)), 217--288, https://epubs.siam.org/doi/abs/10.1137/090771806
"""
raise NotImplementedError

[docs]def randproj_dct_real(trans, a, k, statecomm):
r"""
randproj_dct_real computes a fast random projection of a real :math:m\times n matrix :math:A using a Discrete Cosine Transform (DCT).
The function can be used as a building block in Randomised Numerical Linear Algebra (RNLA) algorithms, such as :math:QR decompositions, Singular Value Decompositions (SVDs), and (approximate) least squares solvers.

.. _f10da-py2-py-doc:

For full information please refer to the NAG Library document for f10da

https://www.nag.com/numeric/nl/nagdoc_28.4/flhtml/f10/f10daf.html

.. _f10da-py2-py-parameters:

**Parameters**
**trans** : str, length 1
Specifies whether the operation pre-multiplies :math:A or post-multiplies :math:A.

:math:\mathrm{trans} = \texttt{'N'}

Random projection is done by post-multiplication, :math:Y = A\Omega.

:math:\mathrm{trans} = \texttt{'T'}

Random projection is done by pre-multiplication, :math:Y = \Omega A.

**a** : float, array-like, shape :math:\left(m, n\right)
The :math:m\times n matrix :math:A.

**k** : int
:math:k, number of columns in the random projection, :math:Y = A\Omega.

**statecomm** : dict, RNG communication object, modified in place
RNG communication structure.

This argument must have been initialized by a prior call to :meth:rand.init_repeat <naginterfaces.library.rand.init_repeat> or :meth:rand.init_nonrepeat <naginterfaces.library.rand.init_nonrepeat>.

**Returns**
**a** : float, ndarray, shape :math:\left(:, :\right)
If :math:\mathrm{trans} = \texttt{'N'}, the :math:m\times k random projection of :math:A.

If :math:\mathrm{trans} = \texttt{'T'}, the :math:k\times n random projection of :math:A.

.. _f10da-py2-py-errors:

**Raises**
**NagValueError**
(errno :math:1)
On entry, :math:\mathrm{trans} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{trans} = \texttt{'N'} or :math:\texttt{'T'}.

(errno :math:2)
On entry, :math:m = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:m > 0.

(errno :math:3)
On entry, :math:n = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:n > 0.

(errno :math:4)
On entry, :math:\mathrm{k} = \langle\mathit{\boldsymbol{value}}\rangle and :math:\mathrm{trans} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: if :math:\mathrm{trans} = \texttt{'N'}, :math:0 < \mathrm{k}\leq n, if :math:\mathrm{trans} = \texttt{'T'}, :math:0 < \mathrm{k}\leq m.

(errno :math:6)
On entry, :math:\mathrm{statecomm}\ ['state'] vector has been corrupted or not initialized.

.. _f10da-py2-py-notes:

**Notes**
A random projection is written as either

.. math::
Y = A\Omega

or

.. math::
Y = \Omega A\text{,}

where :math:A is a real :math:m\times n matrix, :math:\Omega is an :math:n\times k matrix in the first case, and a :math:k\times m matrix in the second case.
These cases are referred to as random projection by post-multiplication and random projection by pre-multiplication, respectively.

When a random projection by post-multiplication uses the DCT, it is written as

.. math::
\Omega = DFC\text{,}

where :math:D is a diagonal matrix whose values are uniformly distributed on the set :math:\left\{-1, +1\right\}, :math:F is a DCT, and :math:C is a matrix that selects a subset of columns, uniformly at random.

When a random projection by pre-multiplication uses the DCT, it is written as

.. math::
\Omega = RFD\text{.}

The operators :math:D and :math:F are as above and :math:R is a matrix that selects a subset of rows, again uniformly at random.

None of these matrix operators require a full matrix-matrix product to be computed.
The computational complexity of applying this type of random projection is :math:\mathrm{O}\left(mn\log\left(k\right)\right).
More details of the DCT-based random projection can be found in Avron et al. (2010).

The DCT-based random projection is closely related to the Subsampled Random Fourier Transform (SRFT) presented in Section 4.6 of Halko et al. (2011).

.. _f10da-py2-py-references:

**References**
Avron, H, Maymounkov, P and Toledo, S, 2010, Blendenpik: Supercharging LAPACK's least-squares solver, SIAM J. Sci. Comput. (32(3)), 1217--1236

Halko, N, 2012, Randomized methods for computing low-rank approximations of matrices, PhD thesis

Halko, N, Martinsson, P G and Tropp, J A, 2011, Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions, SIAM Rev. (53(2)), 217--288, https://epubs.siam.org/doi/abs/10.1137/090771806
"""
raise NotImplementedError