Decompositions and Solve

In matrix calculus the decomposition of a matrix into the matrix product of several matrices with special properties (eg. into an orthogonal and a triangular matrix (QR) or orthogonal matrices and diagonal of singular values (SVD)) are among the most powerful tools to devise numerical algorithms. In the case of tensors of higher degree it is necessary to indicate along which modes the decomposition is supposed to happen, so xerus uses the notation of indexed equations explained in the previous chapter (Indices and Equations).

QR Decompositions

To provide an intuitive approach to decompositions, xerus uses the assignment of multiple tensors with a single operator to denote them. Here (Q, R) = QR(A) reads as “Q and R are defined as the QR-Decomposition of A”, even though we naturally have to provide indices to make this line well defined:

// perform QR decomposition of A and store result in Q and R
(Q(i,r), R(r,j)) = xerus.QR(A(i,j));

# perform QR decomposition of A and store result in Q and R
(Q(i,r), R(r,j)) << xerus.QR(A(i,j))

In these decompositions we can distribute the modes of the original tensor as we please, but to be well defined, Q and R must only share one and exactly one index.

// well defined QR decomposition of a degree 4 tensor
(Q(i,r,k), R(l,r,j)) = xerus.QR(A(i,j,k,l));
// invalid: (Q(i,r,s), R(r,s,j)) = xerus.QR(A(i,j));

# well defined QR decomposition of a degree 4 tensor
(Q(i,r,k), R(l,r,j)) << xerus.QR(A(i,j,k,l))
# invalid: (Q(i,r,s), R(r,s,j)) << xerus.QR(A(i,j))

For convenience xerus defines four variants of the QR decomposition. Assuming the input is of size $m\times n$, $min = \operatorname{min}(m,n)$ and $r$ the rank of the input we have following resulting objects:

Singular Value Decompositions

The Singular Value Decomposition in xerus is called very much like the QR decomposition:

// calculate the SVD of A and store the resulting matrices in U, S and Vt
(U(i,r1), S(r1,r2), Vt(r2,j)) = xerus.SVD(A(i,j));

# calculate the SVD of A and store the resulting matrices in U, S and Vt
(U(i,r1), S(r1,r2), Vt(r2,j)) << xerus.SVD(A(i,j))

In this form it is rank-revealing (so S is of dimensions $r\times r$ instead of $\operatorname{min}(m,n)\times\operatorname{min}(m,n)$) and exact, but it is possible to pass optional arguments to use it as a truncated SVD.

// calculate the SVD, truncated to at most 5 singular values
size_t numSingularVectors = 5;
// or until a singular value is smaller than 0.01 times the maximal singular value
double epsilon = 0.01;
(U(i,r1), S(r1,r2), Vt(r2,j)) = xerus.SVD(A(i,j), numSingularVectors, epsilon);

# calculate the SVD, truncated to at most 5 singular values
numSingularVectors = 5
# or until a singular value is smaller than 0.01 times the maximal singular value
epsilon = 0.01
(U(i,r1), S(r1,r2), Vt(r2,j)) = xerus.SVD(A(i,j), maxRank=numSingularVectors, eps=epsilon);

Solving Linear Equations

A common application for matrix decompositions is to solve matrix equations of the form $A\cdot x = b$ for $x$ via QR, LU, LDL$^T$ or Cholesky decompositions. In xerus this is again provided in indexed notation via the operator/.

// solve A(i,j)*x(j) = b(i) for x
x(j) = b(i) / A(i,j);

# solve A(i,j)*x(j) = b(i) for x
x(j) << b(i) / A(i,j)

Depending on the representation and some properties of A, xerus will automatically choose one of the above mentioned decompositions to solve the system.