# -*- coding: utf-8 -*-
r"""
Module Summary
--------------
Interfaces for the NAG Mark 30.2 `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://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/f10/f10intro.html
"""
# NAG Copyright 2017-2024.
[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://support.nag.com/numeric/nl/nagdoc_30.2/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://support.nag.com/numeric/nl/nagdoc_30.2/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://support.nag.com/numeric/nl/nagdoc_30.2/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
