Approximation

To convert dense data into low-rank format, different approximation algorithms have been developed over the last years in addition to the classical approximation algorithms. HLR implements several of those, with the main focus on approximation during H-matrix arithmetic, i.e., aside from the direct approximation of dense matrices, also (re-) approximation of given low-rank data or general operators is available.

Singular Value Decomposition

For each matrix \(M\) there exists a decomposition

\[ M = U·S·V^H \]

with orthogonal \(U\) and \(V\) and diagonal \(S = diag(s_1,...,s_n)\) with singular values \(s_1 ≥ s_2 ≥ ... ≥ s_n\).

By choosing \(k ≤ n\) the truncation

\[ U(1:k,:)·S(1:k,1:k)·V(1:k,:)^H \]

one obtains the best approximation to \(M\) with respect to all rank-\(k\) matrices.

The SVD can also efficiently be computed for low-rank matrices \(W·X^H\), e.g., to truncate from rank \(k\) to rank \(k' < k\), with the same approximation properties:

function [ W, X ] = svd ( U, V, eps )
  [ QU, RU ] = qr( U, 0 );
  [ QV, RV ] = qr( V, 0 );
  R = RU * RV';
  [ Us, Ss, Vs ] = svd( R );
  k = trunc_rank( Ss );
  W = QU*Us(:,1:k)*S(1:k,1:k);
  X = QV*Vs(:,1:k);
end

Rank Revealing QR

In Rank Revealing QR (RRQR), a QR decomposition of a dense matrix is computed with pivoting such that the order of the Housholder reflections represents their significance on the original matrix:

\[ M·P = Q·R \]

This can be used to select the first \(k\) pairs for the low-rank approximation by examining the norm of the matrices \(R(1:i,1:i)\) for \(i = 1,...,n\). Due to this, RRQR permits error control of the approximation, which also enables variable rank approximation.

function [ W, X ] = rrqr ( M, eps )
  # column pivoted QR
  [ Q, R, P ] = qr( M, 0 );

  for i = 1:n
    S(i) = norm( R(i:n,i:n), "fro" );
    if S(i) < eps * S(1)
      k = i;
      break;
    end
  end

  W = Q(:,1:k);
  for i = 1:ncols
    X(P(i),:) = R'(i,1:k);
  end
end

Randomized SVD

By using randomized variables a column basis of rank \(k\) for the given operator is determined by matrix-vector multiplication. Afterwards the operator is then represented within this basis set, resulting in a rank \(k\) low-rank approximation.

With iterative application of this process, also a variable rank approximation is possible.

Due to only matrix-vector multiplication required for the computation, general operators are possible, e.g., implicit matrix sums.

However, because of the random nature of the process, the final error is only guaranteed within a given probability.

function [ W, X ] = randsvd( M, eps )
  B          = column_basis( M, eps );
  [Q,R]      = qr( M'·B, 0 );
  [Us,Ss,Vs] = svd( R );
  k          = rank( Ss, eps );
  W          = B·Vs(:,1:k)·S(1:k,1:k);
  X          = M'·B·Us(:,1:k);
end

Randomized Low-Rank

Randomized Low-Rank approximation is a simplification of Randomized SVD by directly using the column basis for approximation without an addition SVD.

function [ W, X ] = randlr ( M, eps )
  W := column_basis( M, eps );
  X := M'·W;
end

Adaptive Cross Approximation

In ACA, rows and column of the matrix are chosen and used as rank-1 updates to the low-rank approximation. The criterium for the specific choice of these rows/column differs between the variants of ACA, e.g., ACA+, ACA-Full, etc.

Below, the next column is determined by the maximal entry in the current row.

function [ W, X ] = aca( M, eps )
  c_1 = 1;
  for i = 1,n
    w_i = M(:,c_i) - W · X(c_i,:)';
    r_i = maxidx(w_i);
    w_i = w_i / w_i(r_i);
    x_i = M(r_i,:)' - X · W(r_i,:)';
    W   = [ W, w_i ];
    X   = [ X, x_i ];
    if norm( w_i · x_i', "fro" ) <= eps norm( W · X', "fro" )
      break;
    end
    c_i+1 = maxidx(x_i);
  end
end

Lanczos Bidiagonalization

Iteratively computes bases \(W_k\) and \(X_k\) of the Kyrlov spaces

\[ K(MM^H, w_1) = span( \{ (MM^H)^i w_1 : 0 \le i \le k \} ) \]

and $\( K(M^HM, x_1) = span( \{ (M^HM)^i x_1 : 0 \le i \le k \} ) \)$

with random \(w_1\), \(x_1 = M^H w_1 / |M^H w_1|\) such that

\[ M \approx W_k B_k X_k^H \]

with bidiagonal \(B_k\):

\[\begin{split} B_k = \begin{pmatrix} \alpha_1 & & & \\\\ \beta_2 & \alpha\_2 & & \\\\ & \ddots & \ddots & \\\\ & & \beta_k & \alpha_k \end{pmatrix} \end{split}\]

function [ W, X ] = lanczos ( M, eps )
  w_1 = rand( n );
  w_1 = w_1 / norm(w_1);
  x_1 = M' · w_1;
  α_1 = norm( x_1 );
  x_1 = x_1 / α_1;
  W   = [ w_1 ];
  X   = [ x_1 ];
  for i = 2:n
    w_i = M · x_(i-1) - α_(i-1)·w_(i-1);
    β_i = norm( w_i );
    w_i = w_i / β_i;
    w_i = w_i - Σ_(j<i-1) w_j'·w_i·w_j;
    x_i = M'·w_i - β_i·x_(i-1);
    α_i = norm( x_i );
    x_i = x_i / α_i;
    x_i = x_i - Σ_(j<i-1) x_j'·x_i·x_j;
    W   = [ W, w_i ];
    X   = [ X, x_i ];
    if norm( w_i·x_i', "fro" ) ≤ eps·norm( W·X', "fro" )
      break
    end
  end
  W = W·diag( [ α_1, α_2, … ] ) + diag( [ β_2, β_3, … ], -1 );
end