NAG FL Interface
f10daf (randproj_dct_real)
1
Purpose
f10daf computes a fast random projection of a real $m$ by $n$ matrix $A$ using a Discrete Cosine Transform (DCT). The routine can be used as a building block in Randomised Numerical Linear Algebra (RNLA) algorithms, such as $QR$ decompositions, Singular Value Decompositions (SVDs), and (approximate) least squares solvers.
2
Specification
Fortran Interface
Integer, Intent (In) 
:: 
m, n, lda, k 
Integer, Intent (Inout) 
:: 
state(*), ifail 
Real (Kind=nag_wp), Intent (Inout) 
:: 
a(lda,*) 
Character (1), Intent (In) 
:: 
trans 

C Header Interface
#include <nag.h>
void 
f10daf_ (const char *trans, const Integer *m, const Integer *n, double a[], const Integer *lda, const Integer *k, Integer state[], Integer *ifail, const Charlen length_trans) 

C++ Header Interface
#include <nag.h> extern "C" {
void 
f10daf_ (const char *trans, const Integer &m, const Integer &n, double a[], const Integer &lda, const Integer &k, Integer state[], Integer &ifail, const Charlen length_trans) 
}

The routine may be called by the names f10daf or nagf_rnla_randproj_dct_real.
3
Description
A random projection is written as either
or
where
$A$ is a real
$m$ by
$n$ matrix,
$\Omega $ is an
$n$ by
$k$ matrix in the first case, and a
$k$ by
$m$ matrix in the second case. These cases are referred to as random projection by postmultiplication and random projection by premultiplication, respectively.
When a random projection by postmultiplication uses the DCT, it is written as
where
$D$ is a diagonal matrix whose values are uniformly distributed on the set
$\left\{1,\mathrm{+1}\right\}$,
$F$ is a DCT, and
$C$ is a matrix that selects a subset of columns, uniformly at random.
When a random projection by premultiplication uses the DCT, it is written as
The operators
$D$ and
$F$ are as above and
$R$ is a matrix that selects a subset of rows, again uniformly at random.
None of these matrix operators require a full matrixmatrix product to be computed. The computational complexity of applying this type of random projection is
$\mathit{O}\left(mn\mathrm{log}\left(k\right)\right)$. More details of the DCTbased random projection can be found in
Avron et al. (2010).
The DCTbased random projection is closely related to the Subsampled Random Fourier Transform (SRFT) presented in Section 4.6 of
Halko et al. (2011).
4
References
Avron H, Maymounkov P and Toledo S (2010) Blendenpik: Supercharging LAPACK's leastsquares solver SIAM J. Sci. Comput. 32(3) 1217–1236
Halko N (2012) Randomized methods for computing lowrank 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
5
Arguments

1:
$\mathbf{trans}$ – Character(1)
Input

On entry: specifies whether the operation premultiplies
$A$ or postmultiplies
$A$.
 ${\mathbf{trans}}=\text{'N'}$
 Random projection is done by postmultiplication, $Y=A\Omega $.
 ${\mathbf{trans}}=\text{'T'}$
 Random projection is done by premultiplication, $Y=\Omega A$.
Constraint:
${\mathbf{trans}}=\text{'N'}$ or $\text{'T'}$.

2:
$\mathbf{m}$ – Integer
Input

On entry: $m$, the number of rows of the matrix $A$.
Constraint:
${\mathbf{m}}>0$.

3:
$\mathbf{n}$ – Integer
Input

On entry: $n$, the number of columns of the matrix $A$.
Constraint:
${\mathbf{n}}>0$.

4:
$\mathbf{a}\left({\mathbf{lda}},*\right)$ – Real (Kind=nag_wp) array
Input/Output

Note: the second dimension of the array
a
must be at least
${\mathbf{n}}$.
On entry: the $m$ by $n$ matrix $A$.
On exit: if
${\mathbf{trans}}=\text{'N'}$, the first
$k$ columns of
$A$ are overwritten with the
$m$ by
$k$ random projection of
$A$.
If ${\mathbf{trans}}=\text{'T'}$, the first $k$ rows of $A$ are overwritten with the $k$ by $n$ random projection of $A$.

5:
$\mathbf{lda}$ – Integer
Input

On entry: the first dimension of the array
a as declared in the (sub)program from which
f10daf is called.
Constraint:
${\mathbf{lda}}\ge {\mathbf{m}}$.

6:
$\mathbf{k}$ – Integer
Input

On entry: $k$, number of columns in the random projection, $Y=A\Omega $.
Constraints:
 if ${\mathbf{trans}}=\text{'N'}$, $0<{\mathbf{k}}\le {\mathbf{n}}$;
 if ${\mathbf{trans}}=\text{'T'}$, $0<{\mathbf{k}}\le {\mathbf{m}}$.

7:
$\mathbf{state}\left(*\right)$ – Integer array
Communication Array
Note: the actual argument supplied
must be the array
state supplied to the initialization routines
g05kff or
g05kgf.
On entry: contains information on the selected base generator and its current state.
On exit: contains updated information on the state of the generator.

8:
$\mathbf{ifail}$ – Integer
Input/Output

On entry:
ifail must be set to
$0$,
$1\text{or}1$. If you are unfamiliar with this argument you should refer to
Section 4 in the Introduction to the NAG Library FL Interface for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value
$1\text{or}1$ is recommended. If the output of error messages is undesirable, then the value
$1$ is recommended. Otherwise, if you are not familiar with this argument, the recommended value is
$0$.
When the value $\mathbf{1}\text{or}\mathbf{1}$ is used it is essential to test the value of ifail on exit.
On exit:
${\mathbf{ifail}}={\mathbf{0}}$ unless the routine detects an error or a warning has been flagged (see
Section 6).
6
Error Indicators and Warnings
If on entry
${\mathbf{ifail}}=0$ or
$1$, explanatory error messages are output on the current error message unit (as defined by
x04aaf).
Errors or warnings detected by the routine:
 ${\mathbf{ifail}}=1$

On entry, ${\mathbf{trans}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{trans}}=\text{'N'}$ or $\text{'T'}$.
 ${\mathbf{ifail}}=2$

On entry, ${\mathbf{m}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{m}}>0$.
 ${\mathbf{ifail}}=3$

On entry, ${\mathbf{n}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{n}}>0$.
 ${\mathbf{ifail}}=4$

On entry, ${\mathbf{k}}=\u2329\mathit{\text{value}}\u232a$ and ${\mathbf{trans}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: if ${\mathbf{trans}}=\text{'N'}$, $0<{\mathbf{k}}\le {\mathbf{n}}$, if ${\mathbf{trans}}=\text{'T'}$, $0<{\mathbf{k}}\le {\mathbf{m}}$.
 ${\mathbf{ifail}}=6$

On entry,
state vector has been corrupted or not initialized.
 ${\mathbf{ifail}}=7$

On entry, ${\mathbf{lda}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{lda}}\ge {\mathbf{m}}$.
 ${\mathbf{ifail}}=99$
An unexpected error has been triggered by this routine. Please
contact
NAG.
See
Section 7 in the Introduction to the NAG Library FL Interface for further information.
 ${\mathbf{ifail}}=399$
Your licence key may have expired or may not have been installed correctly.
See
Section 8 in the Introduction to the NAG Library FL Interface for further information.
 ${\mathbf{ifail}}=999$
Dynamic memory allocation failed.
See
Section 9 in the Introduction to the NAG Library FL Interface for further information.
7
Accuracy
The accuracy of algorithms that use
f10daf depend on the extent to which the random projection,
$Y$, captures the range (i.e., column space) of the matrix
$A$. More formally, the following probabilistic error bound holds,
with failure probability at most
$\mathit{O}\left({k}^{1}\right)$, and where
${\sigma}_{k+1}$ is the
$\left(k+1\right)$th singular value and the norm on the lefthand side of the equation is the spectral norm.
The matrix $\left(I{P}_{Y}\right)$ is sometimes referred to as the orthogonal projector onto the complementary subspace, ${\mathrm{range}\left({P}_{Y}\right)}^{\perp}$.
Informally,
$Y$ captures the range of
$A$ well if
is dominated by the first term,
${P}_{Y}A$.
See Section 8–11 of
Halko et al. (2011) for more details.
8
Parallelism and Performance
f10daf is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
f10daf makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.
Please consult the
X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the
Users' Note for your implementation for any additional implementationspecific information.
f10daf uses the same DCT as
c06rff. The time taken by
c06rff is fastest if the only prime factors of the projected dimension are
$2$,
$3$ or
$5$. In order to make the performance of
f10daf less sensitive to the size of the projected dimension, the input,
a, is padded with zeros up to the nearest power of
$2$ or the nearest
$1000$, whichever is smaller.
10
Example
This example applies a random projection using the DCT to the
$6$ by
$6$ matrix
It then computes a
$6$ by
$5$ orthornormal matrix
$Q$ and a probabilistic error bound,
$E$ that quantifies how well the range of
$Q$ approximates the range of
$A$. The matrix
$Q$ is calculated through a
$QR$ factorization of the projection
$Y=A\Omega $. The error bound is based on the following result from Section 4.3 of
Halko et al. (2011). Given a sequence
$\left\{{\omega}^{\left(i\right)}:i=1,2,\dots ,r\right\}$ of standard Gaussian vectors, where
$r$ is some positive integer, then the measure,
$E$, of the error in how well the projection captures the range of
$A$ is given by
where
$\Vert \xb7\Vert $ is the spectral norm when applied to a matrix and the Euclidean norm when applied to a vector.
The exact approximation error $\Vert \left(IQ{Q}^{*}\right)A\Vert $ is also computed for comparison.
This example is constructed so that the true rank of $A$ is $4$: the last $3$ rows are identical so the matrix only has $4$ linearly independent rows. This means that when $k\ge 4$ the accuracy of $Q$ will be close to machine precision.
10.1
Program Text
10.2
Program Data
10.3
Program Results