gpp_linear_algebra_test¶

gpp_linear_algebra_test.hpp¶

Functions for testing gpp_linear_algebra and supporting utilities.

Includes Build* functions that build matrices with various properties. For example, * Well-conditioned SPD matrices * Orthogonal matrices * Ill-conditioned SPD

and so forth. These are in turn used as test inputs for the linear algebra routines.

The linear algebra tests are all called through RunLinearAlgebraTests(); these generally involve tests like checking against hand-verified examples, verifying special matrix properties, testing special cases against general cases, and verifying analytic mathematical bounds on residual/error.

namespace optimal_learning

Macro to allow restrict as a keyword for C++ compilation and CUDA/nvcc compilation. See related entry in gpp_common.hpp for more details.

Variables

constexpr double kProlateDefaultParameter

Builds the “Prolate” matrix. This is a generally ill-conditioned, symmetric matrix. With the default parameter, the condition number is about 10^(n/2); e.g., with n=10, the condition number is 1.84 x 10^6.

This matrix is also Toeplitz (constant along diagonals).

For 0 < alpha <= 0.5, the prolate matrix is SPD. At alpha = 0, it is singular; at alpha = 0.5, it is the identity.

This matches the result of: gallery('prolate', size, alpha) in MATLAB. As in MATLAB, the default parameter is 0.25. See implementation for matrix definition.

Parameters:

alpha:	prolate parameter
size:	dimension of prolate matrix

Outputs:

prolate_matrix[size][size]:
	the prolate matrix with parameter alpha, dimension size

constexpr double kMolerDefaultParameter

Builds the “Moler” matrix. This is a generally ill-conditioned, SPD matrix. Roughly speaking, with the default parameter, the condition number is about 10^(n/2); e.g., with n=10, the condition number is 3.68 x 10^6.

Note

some parameter choices (near 0) result in near-identity matrices.

This matrix has generally has small eigenvalue, with the others clustered in a comparatively small range away from the smallest.

This matches the result of: gallery('moler', size, alpha) in MATLAB. As in MATLAB, the default parameter is -1.0. See implementation for matrix definition.

Parameters:

alpha:	moler matrix parameter
size:	dimension of moler matrix

Outputs:

moler_matrix[size][size]:
	the moler matrix with parameter alpha, dimension size

gpp_linear_algebra_test.cpp¶

Routines to test the functions in gpp_linear_algebra.cpp and gpp_linear_algebra-inl.hpp.

This includes a battery of tests that verify that all of our linear algebra subroutines are working correctly. These tests fall into a few categories:

1. manually verifying a general function (e.g., GeneralMatrixVectorMultiply) and then using that to verify special cases (e.g., Triangular and Symmetric multiply).
2. asserting properties of the underlying matrices or operators: e.g., Q*Q^T = I for Q known to be orthogonal, or X*X^-1 = I, etc.
3. checking correctness against simple, hand-verified cases
4. checking results against analytically known (usually norm-wise) error bounds, taking conditioning into account over well- and ill-conditioned inputs

These tests live in the functions called through RunLinearAlgebraTests().

This file also has implementations various Build.*() functions, which provide interesting inputs for the linear algebra testing. These include various random matrices (random, symmetric, SPD), interesting “standard” matrix examples (prolate, moler, orthogonal symmetric) as well as some matrices with interesting properties resulting from things like Householder Reflections. We also have some utilites for data manipulation (e.g., extracting the lower triangle) as well as routines to manipulate condition number (e.g., adding diagonal dominance).

namespace optimal_learning

Macro to allow restrict as a keyword for C++ compilation and CUDA/nvcc compilation. See related entry in gpp_common.hpp for more details.

Functions

void BuildIdentityMatrix(int size_m, double *restrict matrix)

Utility function to generate a m x m identity matrix. matrix is overwritten.

Parameters:

size_m:	dimension of matrix
matrix[size_m][size_m]:
	pointer to space for identity matrix

Outputs:

matrix[size_m][size_m]:
	identity matrix of order size_m

void BuildProlateMatrix(double alpha, int size, double *restrict prolate_matrix)

A_{i,j} = { 2 * \alpha                                 if i == j `` { sin(2 * pi * alpha * k)/ (pi * k) otherwise`` where k = |i - j|

void BuildMolerMatrix(double alpha, int size, double *restrict moler_matrix)

A_{ij} = { i * alpha^2 + 1.0               if i == j `` { min(i, j) * alpha^2 + alpha otherwise``

void BuildOrthogonalSymmetricMatrix(int size, double *restrict orthog_symm_matrix)

Builds a matrix Q s.t. Q * Q^T = I, Q = Q^T, and \|Q*x\| = \|x\|. In particular, Q` is (real) orthogonal AND symmetric.  This is not the only ``Q with the given properties. Q is not SPD.

This is the eigenvector matrix for a n-point second-difference matrix (e.g., discrete hessian).

Builds an orthogonal (unitary) and symmetric matrix (thus involutary): Q * Q^T = I and Q = Q^T.

Q is NOT SPD.

Parameters:

size:	dimension of matrix

Outputs:

orthog_symm_matrix[size][size]:
	size x size matrix satisfying the above properties

void BuildRandomSymmetricMatrix(int size, double left_bound, double right_bound, UniformRandomGenerator * uniform_generator, double *restrict symmetric_matrix)

Randomly generates half (diagonal and one triangle) of a matrix and copies those values into the other half. The result is not guaranteed to have any special properties (e.g., SPD) beyond symmetry.

Builds a symmetric matrix with random entries. Entries are chosen uniformly at random in the range [left_bound, right_bound].

Matrix condition number is generally ~ 10 * size. HOWEVER since this is a random matrix, extremely poor conditioning is possible (even singular matrices are possible, however unlikely).

Parameters:

size:	dimension of matrix
left_bound:	lower bound for matrix entries
right_bound:	upper bound for matrix_entries
uniform_generator[1]:
	a UniformRandomGenerator object providing the random engine for uniform random numbers

Outputs:

uniform_generator[1]:
	UniformRandomGenerator object will have its state changed due to random draws
symmetric_matrix[size][size]:
	a random symmetric matrix

void BuildRandomLowerTriangularMatrix(int size, UniformRandomGenerator * uniform_generator, double *restrict lower_triangular_matrix)

Randomly generates a lower triangular matrix, zeroing the upper triangle.

Builds a lower triangular matrix with random entries; entries \in [0,1]. The strict upper triangle is set to 0.

Matrix condition number is generally ~ 10 * size. HOWEVER since this is a random matrix, extremely poor conditioning is possible (even singular matrices are possible, however unlikely).

Building an SPD matrix by: A = BuildRandomLowerTriangularMatrix SPD = A * A^T usually leads to ill-conditioned SPD matrices because cond(SPD) is cond(A)^2.

Parameters:

size:	dimension of matrix
uniform_generator[1]:
	a UniformRandomGenerator object providing the random engine for uniform random numbers

Outputs:

uniform_generator[1]:
	UniformRandomGenerator object will have its state changed due to random draws
lower_triangular_matrix[size][size]:
	a random, lower triangular matrix

void BuildRandomSPDMatrix(int size, UniformRandomGenerator * uniform_generator, double *restrict spd_matrix)

A matrix A is SPD if and only if it can be cholesky-factored: L * L^T = A. Generate L randomly and form A.

Builds an SPD matrix with random entries; entries \in [0,1].

Matrix condition number is generally 100 * size^2 with high variance. Since this is a random matrix, extremely poor conditioning is possible (even singular matrices are possible, however unlikely).

Builds SPD matrix by computing L * L^T, where L comes from BuildRandomLowerTriangularMatrix().

Parameters:

size:	dimension of matrix
uniform_generator[1]:
	a UniformRandomGenerator object providing the random engine for uniform random numbers

Outputs:

uniform_generator[1]:
	UniformRandomGenerator object will have its state changed due to random draws
spd_matrix[size][size]:
	a random, SPD matrix

void BuildHouseholderReflectorMatrix(double const *restrict vector, int size, double *restrict householder)

The matrix F (householder) as a function of x (vector) where: F = I - 2 * v*v^T, where v = w / ||w||_2, and w = sign(x[0]) * \|x\|_2 * e_0 + x (e_0 is the cartesian unit vector, [1; zeros(n-1,1)])

Generates matrix F (as a function of x) such that F * x = [ ||x||_2; zeros(n-1,1) ] F is orthogonal (unitary). F is symmetric but never SPD.

See implementation for the definition of F.

The eigenvalues of F are +/- 1; 1 has multiplicity n-1, -1 has multiplicity 1 (hence not SPD).

Parameters:

vector[size]:	vector x used to construct F so that F * x = ||x||_2 * e_0
size:	dimension of matrix, vector

Outputs:

householder[size][size]:
	householder matrix, F

void BuildRandomVector(int size, double left_bound, double right_bound, UniformRandomGenerator * uniform_generator, double *restrict vector)

Builds a vector with random entries. Entries \in [left_bound, right_bound].

A random (non-symmetric) m X n matrix can be built too: BuildRandomVector(m*n, matrix);

Parameters:

size:	number of elements in vector
uniform_generator[1]:
	a UniformRandomGenerator object providing the random engine for uniform random numbers

Outputs:

uniform_generator[1]:
	UniformRandomGenerator object will have its state changed due to random draws
vector[size]:	array filled with random entries

bool CheckMatrixIsSymmetric(double const *restrict matrix, int size, double tolerance)

Checks if the input matrix is symmetric to within tolerance.

If the matrix was explicitly constructed to be symmetric (e.g., A_{j,i} copied from A_{i,j}), then tolerance should be set to 0.0. Otherwise tolerance should usually be near std::numeric_limits<double>::epsilon().

Parameters:

matrix[size][size]:
	matrix to check for symmetry
size:	dimension of matrix
tolerance:	amount of element-wise non-symmetry allowed

Returns:

true if matrix is symmetric

int RunLinearAlgebraTests()

Runs a battery of tests on (supporting) linear algebra routines: * cholesky factorization * Solving A * x = b when A is SPD * Formation of A^-1 * y = A * x, general A * C = A * B, general A, B * y = A * x, A = A^T * y = A * x, A lower triangular * computing A^T`

Returns:

number of test failures: 0 if all is working well.