Imports

tensor, sequninit, decomposition_lapack, triangular, auxiliary_lapack

Procs

proc lu_permuted[T: SupportedDecomposition](a: Tensor[T]): tuple[ PL, U: Tensor[T]]: Compute the pivoted LU decomposition of an input matrix a.

The decomposition solves the equation: A = P L U

where:

P is a permutation matrix

L is a lower-triangular matrix with unit diagonal

U is an upper-triangular matrix

Input:

a, a MxN matrix

Output: with K = min(M, N)

PL, the product of P and L, of shape M, K

U, upper-triangular matrix of shape K, N

Source Edit
proc qr[T: SupportedDecomposition](a: Tensor[T]): tuple[Q, R: Tensor[T]]: Compute the QR decomposition of an input matrix a Decomposition is done through the Householder method without pivoting.

Input:

a, matrix of shape M, N

We note K = min(M, N)

Returns:

Q orthonormal matrix of shape M, K

R upper-triangular matrix of shape K, N

Source Edit
proc svd[T: SupportedDecomposition; U: SomeFloat](A: Tensor[T]; _: typedesc[U]): tuple[ U: Tensor[T], S: Tensor[U], Vh: Tensor[T]]: Compute the Singular Value Decomposition of an input matrix a Decomposition is done through recursive divide & conquer.

Input:

A, matrix of shape M, N

Returns: with K = min(M, N)

U: Unitary matrix of shape M, K with left singular vectors as columns

S: Singular values diagonal of length K in decreasing order

Vh: Unitary matrix of shape K, N with right singular vectors as rows

SVD solves the equation: A = U S V.h

with S being a diagonal matrix of singular values

with V being the right singular vectors and V.h being the hermitian conjugate of V for real matrices, this is equivalent to V.t (transpose)

⚠️: Input must not contain NaN

Compared to Numpy svd procedure, we default to "full_matrices = false".

Exception:

This can throw if the algorithm did not converge.

Source Edit
proc svd[T: SupportedDecomposition](A: Tensor[T]): auto: Source Edit
proc symeig[T: SupportedDecomposition](a: Tensor[T]; return_eigenvectors: static bool = false; uplo: static char = 'U'): tuple[eigenval, eigenvec: Tensor[T]] {. inline.}: Compute the eigenvalues and eigen vectors of a symmetric matrix Input:
A symmetric matrix of shape n x n

A boolean: true if you also want the eigenvectors, false otherwise

A char U for upper or L for lower This allows you to only fill half of the input symmetric matrix

Returns:

A tuple with:
The eigenvalues sorted from lowest to highest. (shape n)

The corresponding eigenvectors of shape n, n if it was requested. If eigenvectors were not requested, this have to be discarded. Using the result will create a runtime error.

Implementation is done through the Multiple Relatively Robust Representations
Source Edit
proc symeig[T: SupportedDecomposition](a: Tensor[T]; return_eigenvectors: static bool = false; uplo: static char = 'U'; slice: HSlice[[type node], [type node]]): tuple[ eigenval, eigenvec: Tensor[T]] {.inline.}: Compute the eigenvalues and eigen vectors of a symmetric matrix Input:
A symmetric matrix of shape n, n

A boolean: true if you also want the eigenvectors, false otherwise

A char U for upper or L for lower This allows you to only fill half of the input symmetric matrix

A slice of the rankings of eigenvalues you request. For example requesting eigenvalues 2 and 3 would be done with 1..2.

Returns:

A tuple with:
The eigenvalues sorted from lowest to highest. (shape m where m is the slice size)

The corresponding eigenvector if it was requested. (shape n, m) If eigenvectors were not requested, this have to be discarded. Using the result will create a runtime error.

Implementation is done through the Multiple Relatively Robust Representations
Source Edit