Matrix operations

This section is concerned with the matrix computations with HODLR matrix.

Note

All the following, we recover the HODLR representation to dense matrix via dense(.), the alternative way is to use recover(.) or simply use the member function of the HODLR representation via H.dense or H.dense().

mhodlr enables multiple precision control for matrix computation based on HODLR representation. It allows users to control the precision in a global environment, without the need to specify the precision everytime when the function is called.

The precision setting is performed by

u = precision('h'); % or ther precision customization
set_prec(u);

Matrix transpose and inverse

The transpose of the HODLR matrix can be simply performed by

tA = hA.transpose();
disp(norm(hadd(tA, A', '-', 'dense'), 'fro')); % print error

The inverse of the HODLR matrix is

rng(0); %fix randomness
A = rand(500);
depth = 5;
min_block_size = 2;
epsilon = 1e-8;
hA = hodlr(A, depth, min_block_size, 'svd', epsilon); % or simply use ``hA = hodlr(A)`` by omitting other parameters as default
iA = inverse(hA);  % return an the inverse in hodlr format
disp(norm(hdot(hA, iA, 'dense') - eye(500), 'fro'))  % print error

The second parameter will be defaulted as ‘hodlr’, which determine the output as hodlr format; One can set it to ‘dense’ to return a dense matrix, for example:

iA = inverse(hA, 'dense');  % return an the inverse in dense format
disp(norm(iA * A - eye(500), 'fro')); % print error

We provide the optional parameters to specify the output format, namely, hodlr and dense, and the third parameter is the algorithm to perform.

Regarding the specified-precision inverse, one can use

u = precision('s');
iA = minverse(hA, u);
disp(norm(hdot(hA, iA, 'dense') - eye(50), 'fro'))

iA = minverse(hA, 'dense');  % return an the inverse in dense format
disp(norm(iA * A - eye(50), 'fro'))

which is similar to the native inverse, but the second parameter is specifed as the precision, while the other two remain the same.

Note

Note that we provide explicit way to truncate HODLR matrix to rank-p HODLR matrix, to do that, one can simply perform

hA = htruncate(hA, 3) # return rank-3 HODLR matrix

Note

The HODLR matrix arithmetic for summation, subtraction, and multiplication can be completed via the shortcut add(A, B), sub(A, B) and dot(A, B) if one just want the HODLR output. If one want a dense output, please following the interfaces as below.

Matrix summation and subtraction

Working precision

One can simply perform the HODLR summation and substraction via the method add and sub. where the inside parameters can be a varying number of inputs, and output is a HODLR matrix.

C1 = add(hA, B, hA);
C2 = A + B + A;
disp(norm(C1.dense - C2, 'fro'));

C1 = sub(hA, B, hA);
C2 = A - B - A;
disp(norm(C1.dense - C2, 'fro'));

Multiple precision

C1 = madd(hA, B, hA);
C2 = A + B + A;
disp(norm(C1.dense - C2, 'fro'));

C1 = msub(hA, B, hA);
C2 = A - B - A;
disp(norm(C1.dense - C2, 'fro'));

Note

The lower level of the summation and subtraction of HODLR matrices (A + B and A - B, A or B are not necessarily of hodlr format) are performed by the method hadd. The are four dominant parameters, i.e., input matrix A, input matrix B, symbol (‘+’ for summation, and ‘-’ for subtraction) denoting summation and subtraction, and output format (‘dense’ or ‘hodlr’, the ‘hodlr’ is default). To perform the operation of summation (subtraction is similar), one can use

hadd(A, B, operation('+'/'-'), 'dense'/'hodlr')

For example:

rng(0); %fix randomness
A = rand(500);
B = rand(500);

depth = 5;
min_block_size = 2;
epsilon = 1e-8;
hA = hodlr(A, depth, min_block_size, 'svd', epsilon); % or simply use ``hA = hodlr(A)`` by omitting other parameters as default
hB = hodlr(B, depth, min_block_size, 'svd', epsilon); % or simply use ``hA = hodlr(A)`` by omitting other parameters as default

disp(norm(hadd(hA, B, '-', 'dense') - (A-B), 'fro'));
disp(norm(hadd(hA, hB, '-', 'dense') - (A-B), 'fro'));
disp(norm(hadd(A, hB, '-', 'dense') - (A-B), 'fro'));
disp(norm(dense(hadd(hA, B, '-', 'hodlr')) - (A-B), 'fro'));
disp(norm(dense(hadd(hA, hB, '-', 'hodlr')) - (A-B), 'fro'));
disp(norm(dense(hadd(A, hB, '-', 'hodlr')) - (A-B), 'fro'));

disp(norm(hadd(hA, B, '+', 'dense') - (A+B), 'fro'));  % print error
disp(norm(hadd(hA, hB, '+', 'dense') - (A+B), 'fro')); % print error
disp(norm(hadd(A, hB, '+', 'dense') - (A+B), 'fro')); % print error
disp(norm(dense(hadd(hA, B, '+', 'hodlr')) - (A+B), 'fro')); % print error
disp(norm(dense(hadd(hA, hB, '+', 'hodlr')) - (A+B), 'fro')); % print error
disp(norm(dense(hadd(A, hB, '+', 'hodlr')) - (A+B), 'fro')); % print error

Note that one should ensure the two inputs are same structure (e.g., same depth) if they both are HODLR format.

Matrix product

Matrix-vector product and matrix-matrix product share the same rountine, one simply use hdot for working precision or mhdot for varying precision to perform comptation.

Working precision

The code example for working precision is as below:

rng(0);
A = rand(100);
x = rand(100, 1); % define vector
X = rand(100, 80); % define matrix

% Usual call for full working precision
hA = hodlr(A, 3, 2, 'svd'); % Use maxmium level of 3 and minimum block size of 2, and perform SVD (default) low rank approximation.
rA = dense(hA);
disp(norm(rA - A, 2)); % print error

b = hdot(hA, x);
err = norm(dense(b) - A * x, 'fro');
disp(err); % print error

b = hdot(hA, x, 'dense');
err = norm(b - A * x, 'fro');
disp(err); % print error

B = hdot(hA, X);
err = norm(dense(B) - A * X, 'fro');
disp(err); % print error

B = hdot(hA, X, 'dense');
err = norm(B - A * X, 'fro');
disp(err); % print error

The third parameter is optional, which indicates whether or not the output is of hodlr format, one can also specify the parameter to dense. The holdr format sometimes requires to be receovered for further operation.

Multiple precision

To simulate specific precision for matrix-matrix product or matrix-vector product, the above example can be simply modifed to:

u = precision('h');

rng(0);
A = rand(100);
x = rand(100, 1); % define vector
X = rand(100, 80); % define matrix

% Usual call for full working precision
hA = hodlr(A, 3, 2, 'svd'); % Use maxmium level of 3 and minimum block size of 2, and perform SVD (default) low rank approximation.

b = mhdot(hA, x, 'dense');
err = norm(b - A * x, 'fro');
disp(err); % print error

B = mhdot(hA, X);
err = norm(dense(B) - A * X, 'fro');
disp(err); % print error

B = mhdot(hA, X, 'dense');
err = norm(B - A * X, 'fro');
disp(err); % print error

LU factorization

Working precision

The LU factorization is performed by the rountine routine

% Output the factors L and U are hodlr format as default
[L, U] = hlu(hA);
err = norm(hdot(L, U, 'dense') - A, 'fro');
disp(err);

% Set the factors L and U to the dense matrix format.
[L, U] = hlu(hA, 'dense');
err = norm(L * U - A, 'fro');
disp(err); % print error

Same to the hdot, the last parameter are used to specify whether or not the output are of hodlr format.

Note

Note that the factors L and U are block lower and upper triangular matrix.

Multiple precision

The working preicion for LU factorization can be specified by the method mhlu:

u = precision('h');
[L, U] = mhlu(hA, 'hodlr');
err = norm(hdot(L, U, 'dense') - A, 'fro');
disp(err); % print error

One can also load the mixed precision mhodlr objects via, for example:

u1 = precision('d');
u2 = precision('s');
u3 = precision('h');
u4 = precision('b');

u_chain = prec_chain(u1, u2, u3, u4);
depth=5;
eps=1e-5;
aphA = amphodlr(u_chain, A, depth, 10, 'svd', eps);
mphA = mphodlr(u_chain, A, depth, 10, 'svd', eps);

u = precision('h'); % set the working precision to half
[L, U] = mhlu(mphA, 'hodlr');
err = norm(hdot(L, U, 'dense') - A, 'fro');
disp(err);

u = precision('s'); % set the working precision to single
[L, U] = mhlu(aphA, 'hodlr');
err = norm(hdot(L, U, 'dense') - A, 'fro');
disp(err); % print error

Cholesky factorization

The Cholesky factorization can be used similarly to LU factorization, which is implemented by the method hchol. The following example briefly illustrate the usage of hchol.

Working precision

rng(0);
R = rand(100);
A = R'*R; % Generate symmetric positive definite matrix

% Usual call for full working precision
hA = hodlr(A, 3, 2, 'svd'); % Use maxmium level of 3 and minimum block size of 2, and perform SVD (default) low rank approximation.

R = hchol(hA); % return a block upper triangular HODLR matrix
disp(norm(hdot(R.transpose(), R, 'dense') - A, 'fro')) % print error

The first and second input of hchol is the input HODLR matrix and format of output, respectively; The second input is optional, which is defaulted as hodlr if it is missing.

To generate the dense output, simply use:

R = hchol(hA, 'dense'); % return a
dusp(norm(R'*R - A, 'fro')); % print error

Multiple precision

The usage of mhchol is similar, it proceeds by simply adding one additional parameter to determine the user-specific working precision.

% Create precisions for each level; Level 1 use precision u1, level 2 use precision u2, ...
u1 = precision('q43');
u2 = precision('q52');
u3 = precision('b');
u4 = precision('s');
u_chain = prec_chain(u1, u2, u3, u4);


% Call mixed precision HODLR representation
amphA = amphodlr(u_chain, A, 3, 2, 'svd'); % Use maxmium level of 3 and minimum block size of 2, and perform SVD (default) low rank approximation.
amprA = dense(amphA);
norm(amprA - A,2) % Compute the error

R = mhchol(amphA); % or R = mhchol(hA);
disp(norm(hdot(R.transpose(), R, 'dense') - A, 'fro')); % print error

Matrix QR factorization

Working precision

We provide three implementations for QR factorizations.

[Q, R] = hqr(hA, 'lt');
[Q, R] = hqr(hA, 'lt2');
[Q, R] = hqr(hA, 'dk');

Multiple precision

[Q, R] = hqr(hA, 'lt');
[Q, R] = hqr(hA, 'lt2');
[Q, R] = hqr(hA, 'dk');

Linear Solver

Solve AX = B

    lu_solve(H, B)
qr_solve(method, H, B) % method is referred to as the method to perform QR factorization