f10daf computes a fast random projection of a real matrix using a Discrete Cosine Transform (DCT). The routine can be used as a building block in Randomised Numerical Linear Algebra (RNLA) algorithms, such as decompositions, Singular Value Decompositions (SVDs), and (approximate) least squares solvers.
The routine may be called by the names f10daf or nagf_rnla_randproj_dct_real.
A random projection is written as either
where is a real matrix, is an matrix in the first case, and a 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
where is a diagonal matrix whose values are uniformly distributed on the set , is a DCT, and 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
The operators and are as above and 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 . 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).
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
On entry: specifies whether the operation pre-multiplies or post-multiplies .
Random projection is done by post-multiplication, .
Random projection is done by pre-multiplication, .
2: – IntegerInput
On entry: , the number of rows of the matrix .
3: – IntegerInput
On entry: , the number of columns of the matrix .
4: – Real (Kind=nag_wp) arrayInput/Output
Note: the second dimension of the array a
must be at least
On entry: the matrix .
On exit: if , the first columns of are overwritten with the random projection of .
If , the first rows of are overwritten with the random projection of .
5: – IntegerInput
On entry: the first dimension of the array a as declared in the (sub)program from which f10daf is called.
6: – IntegerInput
On entry: , number of columns in the random projection, .
if , ;
if , .
7: – Integer arrayCommunication 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: – IntegerInput/Output
On entry: ifail must be set to , or to set behaviour on detection of an error; these values have no effect when no error is detected.
A value of causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of means that an error message is printed while a value of means that it is not.
If halting is not appropriate, the value or is recommended. If message printing is undesirable, then the value is recommended. Otherwise, the value is recommended. When the value or is used it is essential to test the value of ifail on exit.
On exit: unless the routine detects an error or a warning has been flagged (see Section 6).
6Error Indicators and Warnings
If on entry or , explanatory error messages are output on the current error message unit (as defined by x04aaf).
Errors or warnings detected by the routine:
On entry, . Constraint: or .
On entry, . Constraint: .
On entry, . Constraint: .
On entry, and . Constraint: if , , if , .
On entry, state vector has been corrupted or not initialized.
On entry, . Constraint: .
An unexpected error has been triggered by this routine. Please
See Section 7 in the Introduction to the NAG Library FL Interface for further information.
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.
Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.
The accuracy of algorithms that use f10daf depend on the extent to which the random projection, , captures the range (i.e., column space) of the matrix . More formally, the following probabilistic error bound holds,
with failure probability at most , and where is the th singular value and the norm on the left-hand side of the equation is the spectral norm.
The matrix is sometimes referred to as the orthogonal projector onto the complementary subspace, .
Background information to multithreading can be found in the Multithreading documentation.
f10daf is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
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 implementation-specific 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 , or . 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 or the nearest , whichever is smaller.
This example applies a random projection using the DCT to the matrix
It then computes a orthornormal matrix and a probabilistic error bound, that quantifies how well the range of approximates the range of . The matrix is calculated through a factorization of the projection . The error bound is based on the following result from Section 4.3 of Halko et al. (2011). Given a sequence of standard Gaussian vectors, where is some positive integer, then the measure, , of the error in how well the projection captures the range of is given by
where is the spectral norm when applied to a matrix and the Euclidean norm when applied to a vector.
The exact approximation error is also computed for comparison.
This example is constructed so that the true rank of is : the last rows are identical so the matrix only has linearly independent rows. This means that when the accuracy of will be close to machine precision.