Library documentation

abstract type AbstractAlgorithm end

Supertype to dispatch on specific implementations of different the different functions. Concrete subtypes should represent both a way to dispatch to a given implementation, as well as the configuration of that implementation.

See also select_algorithm.

Algorithm{name,KW} <: AbstractAlgorithm

Bare-bones implementation of an algorithm, where name should be a Symbol to dispatch on, and KW is typically a NamedTuple indicating the keyword arguments.

See also @algdef.

LAPACK_Bisection()

Algorithm type to denote the LAPACK driver for computing the eigenvalue decomposition of a Hermitian matrix, or the singular value decomposition of a general matrix using the Bisection algorithm.

LAPACK_DivideAndConquer()

Algorithm type to denote the LAPACK driver for computing the eigenvalue decomposition of a Hermitian matrix, or the singular value decomposition of a general matrix using the Divide and Conquer algorithm.

LAPACK_Expert()

Algorithm type to denote the expert LAPACK driver for computing the Schur or non-Hermitian eigenvalue decomposition of a matrix.

LAPACK_HouseholderLQ(; blocksize, positive = false)

Algorithm type to denote the standard LAPACK algorithm for computing the LQ decomposition of a matrix using Householder reflectors. The specific LAPACK function can be controlled using the keyword arugments, i.e. ?gelqt will be chosen if blocksize > 1 or ?gelqf will be chosen if blocksize == 1. The keyword positive=true can be used to ensure that the diagonal elements of L are non-negative.

LAPACK_HouseholderQR(; blocksize, positive = false, pivoted = false)

Algorithm type to denote the standard LAPACK algorithm for computing the QR decomposition of a matrix using Householder reflectors. The specific LAPACK function can be controlled using the keyword arugments, i.e. ?geqrt will be chosen if blocksize > 1. With blocksize == 1, ?geqrf will be chosen if pivoted == false and ?geqp3 will be chosen if pivoted == true. The keyword positive=true can be used to ensure that the diagonal elements of R are non-negative.

LAPACK_Jacobi()

Algorithm type to denote the LAPACK driver for computing the singular value decomposition of a general matrix using the Jacobi algorithm.

LAPACK_MultipleRelativelyRobustRepresentations()

Algorithm type to denote the LAPACK driver for computing the eigenvalue decomposition of a Hermitian matrix using the Multiple Relatively Robust Representations algorithm.

LAPACK_QRIteration()

Algorithm type to denote the LAPACK driver for computing the eigenvalue decomposition of a Hermitian matrix, or the singular value decomposition of a general matrix using the QR Iteration algorithm.

LAPACK_Simple()

Algorithm type to denote the simple LAPACK driver for computing the Schur or non-Hermitian eigenvalue decomposition of a matrix.

NoTruncation()

Trivial truncation strategy that keeps all values, mostly for testing purposes.

PolarViaSVD(svdalg)

Algorithm for computing the polar decomposition of a matrix A via the singular value decomposition (SVD) of A. The svdalg argument specifies the SVD algorithm to use.

TruncatedAlgorithm(alg::AbstractAlgorithm, trunc::TruncationAlgorithm)

Generic wrapper type for algorithms that consist of first using alg, followed by a truncation through trunc.

TruncationIntersection(trunc::TruncationStrategy, truncs::TruncationStrategy...)

Composition of multiple truncation strategies, keeping values common between them.

TruncationKeepFiltered(filter::Function)

Truncation strategy to keep the values for which filter returns true.

TruncationKeepSorted(howmany::Int, by::Function, rev::Bool)

Truncation strategy to keep the first howmany values when sorted according to by in increasing (decreasing) order if rev is false (true).

abstract type TruncationStrategy end

Supertype to denote different strategies for truncated decompositions that are implemented via post-truncation.

See also truncate!

copy_input(f, A)

Preprocess the input A for a given function, such that it may be handled correctly later. This may include a copy whenever the implementation would destroy the original matrix, or a change of element type to something that is supported.

MatrixAlgebraKit.default_algorithm(f, A; kwargs...)
MatrixAlgebraKit.default_algorithm(f, ::Type{TA}; kwargs...) where {TA}

Select the default algorithm for a given factorization function f and input A. In general, this is called by select_algorithm if no algorithm is specified explicitly. New types should prefer to register their default algorithms in the type domain.

default_pullback_gaugetol(a)

Default tolerance for deciding to warn if incoming adjoints of a pullback rule has components that are not gauge-invariant.

defaulttol(x)

Default tolerance or precision for a given object, e.g. to decide when it can be considerd to be zero or ignored in some other way, or how accurate some quantity needs to be computed.

eig_full(A; kwargs...) -> D, V
eig_full(A, alg::AbstractAlgorithm) -> D, V
eig_full!(A, [DV]; kwargs...) -> D, V
eig_full!(A, [DV], alg::AbstractAlgorithm) -> D, V

Compute the full eigenvalue decomposition of the square matrix A, such that A * V = V * D, where the invertible matrix V contains the eigenvectors and the diagonal matrix D contains the associated eigenvalues.

and therefore will always return complex eigenvalues and eigenvectors. For the real eigenvalue decomposition of symmetric or hermitian operators, see eigh_full.

See also eig_vals(!) and eig_trunc(!).

eig_full(A; kwargs...) -> D, V
eig_full(A, alg::AbstractAlgorithm) -> D, V
eig_full!(A, [DV]; kwargs...) -> D, V
eig_full!(A, [DV], alg::AbstractAlgorithm) -> D, V

Compute the full eigenvalue decomposition of the square matrix A, such that A * V = V * D, where the invertible matrix V contains the eigenvectors and the diagonal matrix D contains the associated eigenvalues.

and therefore will always return complex eigenvalues and eigenvectors. For the real eigenvalue decomposition of symmetric or hermitian operators, see eigh_full.

See also eig_vals(!) and eig_trunc(!).

eig_trunc(A; kwargs...) -> D, V
eig_trunc(A, alg::AbstractAlgorithm) -> D, V
eig_trunc!(A, [DV]; kwargs...) -> D, V
eig_trunc!(A, [DV], alg::AbstractAlgorithm) -> D, V

Compute a partial or truncated eigenvalue decomposition of the matrix A, such that A * V ≈ V * D, where the (possibly rectangular) matrix V contains a subset of eigenvectors and the diagonal matrix D contains the associated eigenvalues, selected according to a truncation strategy.

and therefore will always return complex eigenvalues and eigenvectors. For the real eigenvalue decomposition of symmetric or hermitian operators, see eigh_full.

See also eig_full(!) and eig_vals(!).

eig_trunc(A; kwargs...) -> D, V
eig_trunc(A, alg::AbstractAlgorithm) -> D, V
eig_trunc!(A, [DV]; kwargs...) -> D, V
eig_trunc!(A, [DV], alg::AbstractAlgorithm) -> D, V

Compute a partial or truncated eigenvalue decomposition of the matrix A, such that A * V ≈ V * D, where the (possibly rectangular) matrix V contains a subset of eigenvectors and the diagonal matrix D contains the associated eigenvalues, selected according to a truncation strategy.

and therefore will always return complex eigenvalues and eigenvectors. For the real eigenvalue decomposition of symmetric or hermitian operators, see eigh_full.

See also eig_full(!) and eig_vals(!).

eig_vals(A; kwargs...) -> D
eig_vals(A, alg::AbstractAlgorithm) -> D
eig_vals!(A, [D]; kwargs...) -> D
eig_vals!(A, [D], alg::AbstractAlgorithm) -> D

Compute the list of eigenvalues of A.

and therefore will always return complex eigenvalues and eigenvectors. For the real eigenvalue decomposition of symmetric or hermitian operators, see eigh_full.

See also eig_full(!) and eig_trunc(!).

eig_vals(A; kwargs...) -> D
eig_vals(A, alg::AbstractAlgorithm) -> D
eig_vals!(A, [D]; kwargs...) -> D
eig_vals!(A, [D], alg::AbstractAlgorithm) -> D

Compute the list of eigenvalues of A.

and therefore will always return complex eigenvalues and eigenvectors. For the real eigenvalue decomposition of symmetric or hermitian operators, see eigh_full.

See also eig_full(!) and eig_trunc(!).

eigh_full(A; kwargs...) -> D, V
eigh_full(A, alg::AbstractAlgorithm) -> D, V
eigh_full!(A, [DV]; kwargs...) -> D, V
eigh_full!(A, [DV], alg::AbstractAlgorithm) -> D, V

Compute the full eigenvalue decomposition of the symmetric or hermitian matrix A, such that A * V = V * D, where the unitary matrix V contains the orthogonal eigenvectors and the real diagonal matrix D contains the associated eigenvalues.

and therefore will retain the eltype of the input for the eigenvalues and eigenvectors. For generic eigenvalue decompositions, see eig_full.

See also eigh_vals(!) and eigh_trunc(!).

eigh_full(A; kwargs...) -> D, V
eigh_full(A, alg::AbstractAlgorithm) -> D, V
eigh_full!(A, [DV]; kwargs...) -> D, V
eigh_full!(A, [DV], alg::AbstractAlgorithm) -> D, V

Compute the full eigenvalue decomposition of the symmetric or hermitian matrix A, such that A * V = V * D, where the unitary matrix V contains the orthogonal eigenvectors and the real diagonal matrix D contains the associated eigenvalues.

and therefore will retain the eltype of the input for the eigenvalues and eigenvectors. For generic eigenvalue decompositions, see eig_full.

See also eigh_vals(!) and eigh_trunc(!).

eigh_trunc(A; kwargs...) -> D, V
eigh_trunc(A, alg::AbstractAlgorithm) -> D, V
eigh_trunc!(A, [DV]; kwargs...) -> D, V
eigh_trunc!(A, [DV], alg::AbstractAlgorithm) -> D, V

Compute a partial or truncated eigenvalue decomposition of the symmetric or hermitian matrix A, such that A * V ≈ V * D, where the isometric matrix V contains a subset of the orthogonal eigenvectors and the real diagonal matrix D contains the associated eigenvalues, selected according to a truncation strategy.

and therefore will retain the eltype of the input for the eigenvalues and eigenvectors. For generic eigenvalue decompositions, see eig_full.

See also eigh_full(!) and eigh_vals(!).

eigh_trunc(A; kwargs...) -> D, V
eigh_trunc(A, alg::AbstractAlgorithm) -> D, V
eigh_trunc!(A, [DV]; kwargs...) -> D, V
eigh_trunc!(A, [DV], alg::AbstractAlgorithm) -> D, V

Compute a partial or truncated eigenvalue decomposition of the symmetric or hermitian matrix A, such that A * V ≈ V * D, where the isometric matrix V contains a subset of the orthogonal eigenvectors and the real diagonal matrix D contains the associated eigenvalues, selected according to a truncation strategy.

and therefore will retain the eltype of the input for the eigenvalues and eigenvectors. For generic eigenvalue decompositions, see eig_full.

See also eigh_full(!) and eigh_vals(!).

eigh_vals(A; kwargs...) -> D
eigh_vals(A, alg::AbstractAlgorithm) -> D
eigh_vals!(A, [D]; kwargs...) -> D
eigh_vals!(A, [D], alg::AbstractAlgorithm) -> D

Compute the list of (real) eigenvalues of the symmetric or hermitian matrix A.

and therefore will retain the eltype of the input for the eigenvalues and eigenvectors. For generic eigenvalue decompositions, see eig_full.

See also eigh_full(!) and eigh_trunc(!).

eigh_vals(A; kwargs...) -> D
eigh_vals(A, alg::AbstractAlgorithm) -> D
eigh_vals!(A, [D]; kwargs...) -> D
eigh_vals!(A, [D], alg::AbstractAlgorithm) -> D

Compute the list of (real) eigenvalues of the symmetric or hermitian matrix A.

and therefore will retain the eltype of the input for the eigenvalues and eigenvectors. For generic eigenvalue decompositions, see eig_full.

See also eigh_full(!) and eigh_trunc(!).

MatrixAlgebraKit.findtruncated(values::AbstractVector, strategy::TruncationStrategy)

Generic interface for finding truncated values of the spectrum of a decomposition based on the strategy. The output should be a collection of indices specifying which values to keep. MatrixAlgebraKit.findtruncated is used inside of the default implementation of truncate! to perform the truncation. It does not assume that the values are sorted. For a version that assumes the values are reverse sorted (which is the standard case for SVD) see MatrixAlgebraKit.findtruncated_sorted.

MatrixAlgebraKit.findtruncated_sorted(values::AbstractVector, strategy::TruncationStrategy)

Like MatrixAlgebraKit.findtruncated but assumes that the values are sorted in reverse order. They are assumed to be sorted in a way that is consistent with the truncation strategy, which generally means they are sorted by absolute value but some truncation strategies allow customizing that. However, note that this assumption is not checked, so passing values that are not sorted in the correct way can silently give unexpected results. This is used in the default implementation of svd_trunc!.

initialize_output(f, A, alg)

Whenever possible, allocate the destination for applying a given algorithm in-place. If this is not possible, for example when the output size is not known a priori or immutable, this function may return nothing.

inv_regularized(a::Number, tol=defaulttol(a))
inv_regularized(A::Matrix, tol=defaulttol(A); isposdef = false, kwargs...)

Compute a smooth regularised inverse (L2 Tikhonov regularisation) of a number or square matrix a.

• For numbers, this is given by inv(hypot(a, tol)).

• For matrices, this is computed using the singular value decomposition and aplying inv_regularized to the singular values. If isposdef = true, the singular value decomposition is equivalent to the (Hermitian) eigenvalue decomposition of A and the latter is used instead.

function inv_safe(a::Number, tol=defaulttol(a))

Compute the inverse of a number a, but return zero if a is smaller than tol.

is_left_isometry(A; isapprox_kwargs...) -> Bool

Test whether a linear map is a left isometry, i.e. A' * A ≈ I. The isapprox_kwargs can be used to control the tolerances of the equality.

See also isisometry and is_right_isometry.

is_right_isometry(A; isapprox_kwargs...) -> Bool

Test whether a linear map is a right isometry, i.e. A * A' ≈ I. The isapprox_kwargs can be used to control the tolerances of the equality.

See also isisometry and is_left_isometry.

isisometry(A; side=:left, isapprox_kwargs...) -> Bool

Test whether a linear map is an isometry, where the type of isometry is controlled by kind:

• side = :left : A' * A ≈ I.
• side = :right : A * A' ≈ I.

The isapprox_kwargs are passed on to isapprox to control the tolerances.

New specializations should overload is_left_isometry and is_right_isometry.

See also isunitary.

isunitary(A; isapprox_kwargs...)

Test whether a linear map is unitary, i.e. A * A' ≈ I ≈ A' * A. The isapprox_kwargs are passed on to isapprox to control the tolerances.

See also isisometry.

iszerotangent(x)

Return true if x is of a type that the different AD engines use to communicate a (co)tangent that is identically zero. By overloading this method, and writing pullback definitions in term of it, we will be able to hook into different AD ecosystems

left_null(A; [kind::Symbol, trunc, alg_qr, alg_svd]) -> N
left_null!(A, [N]; [kind::Symbol, alg_qr, alg_svd]) -> N

Compute an orthonormal basis N for the cokernel of the matrix A of size (m, n), i.e. the nullspace of adjoint(A), such that adjoint(A)*N ≈ 0 and N'*N ≈ I. The keyword argument kind can be used to specify the specific orthogonal decomposition that should be used to factor A, whereas trunc can be used to control the the rank of A via its singular values.

trunc can either be a truncation strategy object or a NamedTuple with fields atol, rtol, and maxnullity.

This is a high-level wrapper and will use one of the decompositions qr! or svd! to compute the orthogonal basis N, as controlled by the keyword arguments.

When kind is provided, its possible values are

• kind == :qr: N is computed using the QR decomposition. This requires isnothing(trunc) and left_null!(A, [N], kind=:qr) is equivalent to qr_null!(A, [N], alg) with a default value alg = select_algorithm(qr_compact!, A; positive=true)

• kind == :svd: N is computed using the singular value decomposition and will contain the left singular vectors corresponding to either the zero singular values if trunc isn't specified or the singular values specified by trunc.

When kind is not provided, the default value is :qr when isnothing(trunc) and :svd otherwise. Finally, finer control is obtained by providing an explicit algorithm using the alg_qr and alg_svd keyword arguments, which will only be used by the corresponding factorization backend. If alg_qr or alg_svd are NamedTuples, a default algorithm is chosen with select_algorithm and the NamedTuple is passed as keyword arguments to that algorithm. alg_qr defaults to (; positive=true) so that by default a positive QR decomposition will be used.

See also right_null(!), left_orth(!), right_orth(!)

left_orth(A; [kind::Symbol, trunc, alg_qr, alg_polar, alg_svd]) -> V, C
left_orth!(A, [VC]; [kind::Symbol, trunc, alg_qr, alg_polar, alg_svd]) -> V, C

Compute an orthonormal basis V for the image of the matrix A of size (m, n), as well as a matrix C (the corestriction) such that A factors as A = V * C. The keyword argument kind can be used to specify the specific orthogonal decomposition that should be used to factor A, whereas trunc can be used to control the precision in determining the rank of A via its singular values.

trunc can either be a truncation strategy object or a NamedTuple with fields atol, rtol, and maxrank.

This is a high-level wrapper and will use one of the decompositions qr_compact!, svd_compact!/svd_trunc!, andleft_polar! to compute the orthogonal basis V, as controlled by the keyword arguments.

When kind is provided, its possible values are

• kind == :qr: V and C are computed using the QR decomposition. This requires isnothing(trunc) and left_orth!(A, [VC]) is equivalent to qr_compact!(A, [VC], alg) with a default value alg = select_algorithm(qr_compact!, A; positive=true)

• kind == :polar: V and C are computed using the polar decomposition, This requires isnothing(trunc) and left_orth!(A, [VC]) is equivalent to left_polar!(A, [VC], alg) with a default value alg = select_algorithm(left_polar!, A)

• kind == :svd: V and C are computed using the singular value decomposition svd_trunc! when a truncation strategy is specified using the trunc keyword argument, and using svd_compact! otherwise. V will contain the left singular vectors and C is computed as the product of the singular values and the right singular vectors, i.e. with U, S, Vᴴ = svd(A), we have V = U and C = S * Vᴴ.

When kind is not provided, the default value is :qr when isnothing(trunc) and :svd otherwise. Finally, finer control is obtained by providing an explicit algorithm for backend factorizations through the alg_qr, alg_polar, and alg_svd keyword arguments, which will only be used if the corresponding factorization is called based on the other inputs. If NamedTuples are passed as alg_qr, alg_polar, or alg_svd, a default algorithm is chosen with select_algorithm and the NamedTuple is passed as keyword arguments to that algorithm. alg_qr defaults to (; positive=true) so that by default a positive QR decomposition will be used.

See also right_orth(!), left_null(!), right_null(!)

left_polar(A; kwargs...) -> W, P
left_polar(A, alg::AbstractAlgorithm) -> W, P
left_polar!(A, [WP]; kwargs...) -> W, P
left_polar!(A, [WP], alg::AbstractAlgorithm) -> W, P

Compute the full polar decomposition of the rectangular matrix A of size (m, n) with m >= n, such that A = W * P. Here, W is an isometric matrix (orthonormal columns) of size (m, n), whereas P is a positive (semi)definite matrix of size (n, n).

See also right_polar(!).

left_polar(A; kwargs...) -> W, P
left_polar(A, alg::AbstractAlgorithm) -> W, P
left_polar!(A, [WP]; kwargs...) -> W, P
left_polar!(A, [WP], alg::AbstractAlgorithm) -> W, P

Compute the full polar decomposition of the rectangular matrix A of size (m, n) with m >= n, such that A = W * P. Here, W is an isometric matrix (orthonormal columns) of size (m, n), whereas P is a positive (semi)definite matrix of size (n, n).

See also right_polar(!).

left_polar_pullback!(ΔA, (W, P), (ΔW, ΔP))

Adds the pullback from the left polar decomposition of A to ΔA given the output (W, P) and cotangent (ΔW, ΔP) of left_polar(A).

lq_compact(A; kwargs...) -> L, Q
lq_compact(A, alg::AbstractAlgorithm) -> L, Q
lq_compact!(A, [LQ]; kwargs...) -> L, Q
lq_compact!(A, [LQ], alg::AbstractAlgorithm) -> L, Q

Compute the compact LQ decomposition of the rectangular matrix A of size (m,n), such that A = L * Q where the matrix Q of size (min(m,n), n) has orthogonal rows spanning the image of A', and the matrix L of size (m, min(m,n)) is lower triangular.

See also lq_full(!).

lq_compact(A; kwargs...) -> L, Q
lq_compact(A, alg::AbstractAlgorithm) -> L, Q
lq_compact!(A, [LQ]; kwargs...) -> L, Q
lq_compact!(A, [LQ], alg::AbstractAlgorithm) -> L, Q

Compute the compact LQ decomposition of the rectangular matrix A of size (m,n), such that A = L * Q where the matrix Q of size (min(m,n), n) has orthogonal rows spanning the image of A', and the matrix L of size (m, min(m,n)) is lower triangular.

See also lq_full(!).

lq_compact_pullback!(ΔA, (L, Q), (ΔL, ΔQ);
                        tol::Real=default_pullback_gaugetol(R),
                        rank_atol::Real=tol,
                        gauge_atol::Real=tol)

Adds the pullback from the LQ decomposition of A to ΔA given the output (L, Q) and cotangent (ΔL, ΔQ) of lq_compact(A; positive=true) or lq_full(A; positive=true).

In the case where the rank r of the original matrix A ≈ L * Q (as determined by rank_atol) is less then the minimum of the number of rows and columns , the cotangents ΔL and ΔQ, only the first r columns of L and the first r rows of Q are well-defined, and also the adjoint variables ΔL and ΔQ should have nonzero values only in the first r columns and rows respectively. If nonzero values in the remaining columns or rows exceed gauge_atol, a warning will be printed.

lq_full(A; kwargs...) -> L, Q
lq_full(A, alg::AbstractAlgorithm) -> L, Q
lq_full!(A, [LQ]; kwargs...) -> L, Q
lq_full!(A, [LQ], alg::AbstractAlgorithm) -> L, Q

Compute the full LQ decomposition of the rectangular matrix A, such that A = L * Q where Q is a square unitary matrix with the same number of rows as A and L is a lower triangular matrix with the same size as A.

See also lq_compact(!).

lq_full(A; kwargs...) -> L, Q
lq_full(A, alg::AbstractAlgorithm) -> L, Q
lq_full!(A, [LQ]; kwargs...) -> L, Q
lq_full!(A, [LQ], alg::AbstractAlgorithm) -> L, Q

Compute the full LQ decomposition of the rectangular matrix A, such that A = L * Q where Q is a square unitary matrix with the same number of rows as A and L is a lower triangular matrix with the same size as A.

See also lq_compact(!).

lq_null(A; kwargs...) -> Nᴴ
lq_null(A, alg::AbstractAlgorithm) -> Nᴴ
lq_null!(A, [Nᴴ]; kwargs...) -> Nᴴ
lq_null!(A, [Nᴴ], alg::AbstractAlgorithm) -> Nᴴ

For a (m, n) matrix A, compute the matrix Nᴴ corresponding the final n - min(m, n) rows oft the unitary Q factor in the full LQ decomposition of A, i.e. the rows that are not present in the Q factor of the compact LQ decomposition. The matrix Nᴴ is such that the isometric matrix N = adjoint(Nᴴ) contains an orthonormal basis for the kernel (null space) of A as its columns, i.e. A * N = 0 or thus A * adjoint(Nᴴ) = 0.

See also qr_full(!) and qr_compact(!).

lq_null(A; kwargs...) -> Nᴴ
lq_null(A, alg::AbstractAlgorithm) -> Nᴴ
lq_null!(A, [Nᴴ]; kwargs...) -> Nᴴ
lq_null!(A, [Nᴴ], alg::AbstractAlgorithm) -> Nᴴ

For a (m, n) matrix A, compute the matrix Nᴴ corresponding the final n - min(m, n) rows oft the unitary Q factor in the full LQ decomposition of A, i.e. the rows that are not present in the Q factor of the compact LQ decomposition. The matrix Nᴴ is such that the isometric matrix N = adjoint(Nᴴ) contains an orthonormal basis for the kernel (null space) of A as its columns, i.e. A * N = 0 or thus A * adjoint(Nᴴ) = 0.

See also qr_full(!) and qr_compact(!).

qr_compact(A; kwargs...) -> Q, R
qr_compact(A, alg::AbstractAlgorithm) -> Q, R
qr_compact!(A, [QR]; kwargs...) -> Q, R
qr_compact!(A, [QR], alg::AbstractAlgorithm) -> Q, R

Compute the compact QR decomposition of the rectangular matrix A of size (m,n), such that A = Q * R where the isometric matrix Q of size (m, min(m,n)) has orthogonal columns spanning the image of A, and the matrix R of size (min(m,n), n) is upper triangular.

See also qr_full(!).

qr_compact(A; kwargs...) -> Q, R
qr_compact(A, alg::AbstractAlgorithm) -> Q, R
qr_compact!(A, [QR]; kwargs...) -> Q, R
qr_compact!(A, [QR], alg::AbstractAlgorithm) -> Q, R

Compute the compact QR decomposition of the rectangular matrix A of size (m,n), such that A = Q * R where the isometric matrix Q of size (m, min(m,n)) has orthogonal columns spanning the image of A, and the matrix R of size (min(m,n), n) is upper triangular.

See also qr_full(!).

qr_compact_pullback!(ΔA, (Q, R), (ΔQ, ΔR);
                        tol::Real=default_pullback_gaugetol(R),
                        rank_atol::Real=tol,
                        gauge_atol::Real=tol)

Adds the pullback from the QR decomposition of A to ΔA given the output (Q,R) and cotangent (ΔQ, ΔR) of qr_compact(A; positive=true) or qr_full(A; positive=true).

In the case where the rank r of the original matrix A ≈ Q * R (as determined by rank_atol) is less then the minimum of the number of rows and columns, the cotangents ΔQ and ΔR, only the first r columns of Q and the first r rows of R are well-defined, and also the adjoint variables ΔQ and ΔR should have nonzero values only in the first r columns and rows respectively. If nonzero values in the remaining columns or rows exceed gauge_atol, a warning will be printed.

qr_full(A; kwargs...) -> Q, R
qr_full(A, alg::AbstractAlgorithm) -> Q, R
qr_full!(A, [QR]; kwargs...) -> Q, R
qr_full!(A, [QR], alg::AbstractAlgorithm) -> Q, R

Compute the full QR decomposition of the rectangular matrix A, such that A = Q * R where Q is a square unitary matrix with the same number of rows as A and R is an upper triangular matrix with the same size as A.

See also qr_compact(!).

qr_full(A; kwargs...) -> Q, R
qr_full(A, alg::AbstractAlgorithm) -> Q, R
qr_full!(A, [QR]; kwargs...) -> Q, R
qr_full!(A, [QR], alg::AbstractAlgorithm) -> Q, R

Compute the full QR decomposition of the rectangular matrix A, such that A = Q * R where Q is a square unitary matrix with the same number of rows as A and R is an upper triangular matrix with the same size as A.

See also qr_compact(!).

qr_null(A; kwargs...) -> N
qr_null(A, alg::AbstractAlgorithm) -> N
qr_null!(A, [N]; kwargs...) -> N
qr_null!(A, [N], alg::AbstractAlgorithm) -> N

For a (m, n) matrix A, compute the matrix N corresponding the final m - min(m, n) columns of the unitary Q factor in the full QR decomposition of A, i.e. the columns that are not present in the Q factor of the compact QR decomposition. The isometric matrix N contains an orthonormal basis for the cokernel of A as its columns, i.e. adjoint(A) * N = 0.

See also lq_full(!) and lq_compact(!).

qr_null(A; kwargs...) -> N
qr_null(A, alg::AbstractAlgorithm) -> N
qr_null!(A, [N]; kwargs...) -> N
qr_null!(A, [N], alg::AbstractAlgorithm) -> N

For a (m, n) matrix A, compute the matrix N corresponding the final m - min(m, n) columns of the unitary Q factor in the full QR decomposition of A, i.e. the columns that are not present in the Q factor of the compact QR decomposition. The isometric matrix N contains an orthonormal basis for the cokernel of A as its columns, i.e. adjoint(A) * N = 0.

See also lq_full(!) and lq_compact(!).

right_null(A; [kind::Symbol, alg_lq, alg_svd]) -> Nᴴ
right_null!(A, [Nᴴ]; [kind::Symbol, alg_lq, alg_svd]) -> Nᴴ

Compute an orthonormal basis N = adjoint(Nᴴ) for the kernel or nullspace of the matrix A of size (m, n), such that A*adjoint(Nᴴ) ≈ 0 and Nᴴ*adjoint(Nᴴ) ≈ I. The keyword argument kind can be used to specify the specific orthogonal decomposition that should be used to factor A, whereas trunc can be used to control the the rank of A via its singular values.

trunc can either be a truncation strategy object or a NamedTuple with fields atol, rtol, and maxnullity.

This is a high-level wrapper and will use one of the decompositions lq! or svd! to compute the orthogonal basis Nᴴ, as controlled by the keyword arguments.

When kind is provided, its possible values are

• kind == :lq: Nᴴ is computed using the (nonpositive) LQ decomposition. This requires isnothing(trunc) and right_null!(A, [Nᴴ], kind=:lq) is equivalent to lq_null!(A, [Nᴴ], alg) with a default value alg = select_algorithm(lq_compact!, A; positive=true)

• kind == :svd: N is computed using the singular value decomposition and will contain the left singular vectors corresponding to the singular values that are smaller than max(atol, rtol * σ₁), where σ₁ is the largest singular value of A.

When kind is not provided, the default value is :lq when isnothing(trunc) and :svd otherwise. Finally, finer control is obtained by providing an explicit algorithm using the alg_lq and alg_svd keyword arguments, which will only be used by the corresponding factorization backend. If alg_lq or alg_svd are NamedTuples, a default algorithm is chosen with select_algorithm and the NamedTuple is passed as keyword arguments to that algorithm. alg_lq defaults to (; positive=true) so that by default a positive LQ decomposition will be used.

See also left_null(!), left_orth(!), right_orth(!)

right_orth(A; [kind::Symbol, trunc, alg_lq, alg_polar, alg_svd]) -> C, Vᴴ
right_orth!(A, [CVᴴ]; [kind::Symbol, trunc, alg_lq, alg_polar, alg_svd]) -> C, Vᴴ

Compute an orthonormal basis V = adjoint(Vᴴ) for the coimage of the matrix A, i.e. for the image of adjoint(A), as well as a matrix C such that A = C * Vᴴ. The keyword argument kind can be used to specify the specific orthogonal decomposition that should be used to factor A, whereas trunc can be used to control the precision in determining the rank of A via its singular values.

trunc can either be a truncation strategy object or a NamedTuple with fields atol, rtol, and maxrank.

This is a high-level wrapper and will use one of the decompositions lq_compact!, svd_compact!/svd_trunc!, and right_polar! to compute the orthogonal basis V, as controlled by the keyword arguments.

When kind is provided, its possible values are

• kind == :lq: C and Vᴴ are computed using the QR decomposition, This requires isnothing(trunc) and right_orth!(A, [CVᴴ]) is equivalent to lq_compact!(A, [CVᴴ], alg) with a default value alg = select_algorithm(lq_compact!, A; positive=true)

• kind == :polar: C and Vᴴ are computed using the polar decomposition, This requires isnothing(trunc) and right_orth!(A, [CVᴴ]) is equivalent to right_polar!(A, [CVᴴ], alg) with a default value alg = select_algorithm(right_polar!, A)

• kind == :svd: C and Vᴴ are computed using the singular value decomposition svd_trunc! when a truncation strategy is specified using the trunc keyword argument, and using svd_compact! otherwise. V = adjoint(Vᴴ) will contain the right singular vectors corresponding to the singular values and C is computed as the product of the singular values and the right singular vectors, i.e. with U, S, Vᴴ = svd(A), we have C = rmul!(U, S) and Vᴴ = Vᴴ.

When kind is not provided, the default value is :lq when isnothing(trunc) and :svd otherwise. Finally, finer control is obtained by providing an explicit algorithm for backend factorizations through the alg_lq, alg_polar, and alg_svd keyword arguments, which will only be used if the corresponding factorization is called based on the other inputs. If alg_lq, alg_polar, or alg_svd are NamedTuples, a default algorithm is chosen with select_algorithm and the NamedTuple is passed as keyword arguments to that algorithm. alg_lq defaults to (; positive=true) so that by default a positive LQ decomposition will be used.

See also left_orth(!), left_null(!), right_null(!)

right_polar(A; kwargs...) -> P, Wᴴ
right_polar(A, alg::AbstractAlgorithm) -> P, Wᴴ
right_polar!(A, [PWᴴ]; kwargs...) -> P, Wᴴ
right_polar!(A, [PWᴴ], alg::AbstractAlgorithm) -> P, Wᴴ

Compute the full polar decomposition of the rectangular matrix A of size (m, n) with n >= m, such that A = P * Wᴴ. Here, P is a positive (semi)definite matrix of size (m, m), whereas Wᴴ is a matrix with orthonormal rows (its adjoint is isometric) of size (n, m).

See also left_polar(!).

right_polar(A; kwargs...) -> P, Wᴴ
right_polar(A, alg::AbstractAlgorithm) -> P, Wᴴ
right_polar!(A, [PWᴴ]; kwargs...) -> P, Wᴴ
right_polar!(A, [PWᴴ], alg::AbstractAlgorithm) -> P, Wᴴ

Compute the full polar decomposition of the rectangular matrix A of size (m, n) with n >= m, such that A = P * Wᴴ. Here, P is a positive (semi)definite matrix of size (m, m), whereas Wᴴ is a matrix with orthonormal rows (its adjoint is isometric) of size (n, m).

See also left_polar(!).

right_polar_pullback!(ΔA, (P, Wᴴ), (ΔP, ΔWᴴ))

Adds the pullback from the left polar decomposition of A to ΔA given the output (P, Wᴴ) and cotangent (ΔP, ΔWᴴ) of right_polar(A).

schur_full(A; kwargs...) -> T, Z, vals
schur_full(A, alg::AbstractAlgorithm) -> T, Z, vals
schur_full!(A, [TZv]; kwargs...) -> T, Z, vals
schur_full!(A, [TZv], alg::AbstractAlgorithm) -> T, Z, vals

Compute the full Schur decomposition of the square matrix A, such that A * Z = Z * T, where the orthogonal or unitary matrix Z contains the Schur vectors and the square matrix T is upper triangular (in the complex case) or quasi-upper triangular (in the real case). The list vals contains the (complex-valued) eigenvalues of A, as extracted from the (quasi-)diagonal of T.

schur_full(A; kwargs...) -> T, Z, vals
schur_full(A, alg::AbstractAlgorithm) -> T, Z, vals
schur_full!(A, [TZv]; kwargs...) -> T, Z, vals
schur_full!(A, [TZv], alg::AbstractAlgorithm) -> T, Z, vals

Compute the full Schur decomposition of the square matrix A, such that A * Z = Z * T, where the orthogonal or unitary matrix Z contains the Schur vectors and the square matrix T is upper triangular (in the complex case) or quasi-upper triangular (in the real case). The list vals contains the (complex-valued) eigenvalues of A, as extracted from the (quasi-)diagonal of T.

schur_vals(A; kwargs...) -> vals
schur_vals(A, alg::AbstractAlgorithm) -> vals
schur_vals!(A, [vals]; kwargs...) -> vals
schur_vals!(A, [vals], alg::AbstractAlgorithm) -> vals

Compute the list of eigenvalues of A by computing the Schur decomposition of A.

See also eig_full(!) and eig_trunc(!).

schur_vals(A; kwargs...) -> vals
schur_vals(A, alg::AbstractAlgorithm) -> vals
schur_vals!(A, [vals]; kwargs...) -> vals
schur_vals!(A, [vals], alg::AbstractAlgorithm) -> vals

Compute the list of eigenvalues of A by computing the Schur decomposition of A.

See also eig_full(!) and eig_trunc(!).

MatrixAlgebraKit.select_algorithm(f, A, alg::AbstractAlgorithm)
MatrixAlgebraKit.select_algorithm(f, A, alg::Symbol; kwargs...)
MatrixAlgebraKit.select_algorithm(f, A, alg::Type; kwargs...)
MatrixAlgebraKit.select_algorithm(f, A; kwargs...)
MatrixAlgebraKit.select_algorithm(f, A, (; kwargs...))

Decide on an algorithm to use for implementing the function f on inputs of type A. This can be obtained both for values A or types A.

If alg is an AbstractAlgorithm instance, it will be returned as-is.

If alg is a Symbol or a Type of algorithm, the return value is obtained by calling the corresponding algorithm constructor; keyword arguments in kwargs are passed along to this constructor.

If alg is not specified (or nothing), an algorithm will be selected automatically with MatrixAlgebraKit.default_algorithm and the keyword arguments in kwargs will be passed to the algorithm constructor. Finally, the same behavior is obtained when the keyword arguments are passed as the third positional argument in the form of a NamedTuple.

sign_safe(s::Number)

Compute the sign of a number s, but return +1 if s is zero so that the result is always a number with modulus 1, i.e. an element of the unitary group U(1).

svd_compact(A; kwargs...) -> U, S, Vᴴ
svd_compact(A, alg::AbstractAlgorithm) -> U, S, Vᴴ
svd_compact!(A, [USVᴴ]; kwargs...) -> U, S, Vᴴ
svd_compact!(A, [USVᴴ], alg::AbstractAlgorithm) -> U, S, Vᴴ

Compute the compact singular value decomposition (SVD) of the rectangular matrix A of size (m, n), such that A = U * S * Vᴴ. Here, U is an isometric matrix (orthonormal columns) of size (m, k), whereas Vᴴ is a matrix of size (k, n) with orthonormal rows and S is a square diagonal matrix of size (k, k), with k = min(m, n).

See also svd_full(!), svd_vals(!) and svd_trunc(!).

svd_compact(A; kwargs...) -> U, S, Vᴴ
svd_compact(A, alg::AbstractAlgorithm) -> U, S, Vᴴ
svd_compact!(A, [USVᴴ]; kwargs...) -> U, S, Vᴴ
svd_compact!(A, [USVᴴ], alg::AbstractAlgorithm) -> U, S, Vᴴ

Compute the compact singular value decomposition (SVD) of the rectangular matrix A of size (m, n), such that A = U * S * Vᴴ. Here, U is an isometric matrix (orthonormal columns) of size (m, k), whereas Vᴴ is a matrix of size (k, n) with orthonormal rows and S is a square diagonal matrix of size (k, k), with k = min(m, n).

See also svd_full(!), svd_vals(!) and svd_trunc(!).

svd_compact_pullback!(ΔA, USVᴴ, ΔUSVᴴ;
                        tol::Real=default_pullback_gaugetol(S),
                        rank_atol::Real = tol,
                        degeneracy_atol::Real = tol,
                        gauge_atol::Real = tol)

Adds the pullback from the SVD of A to ΔA given the output USVᴴ of svd_compact or svd_full and the cotangent ΔUSVᴴ of svd_compact, svd_full or svd_trunc.

In particular, it is assumed that A ≈ U * S * Vᴴ, or thus, that no singular values with magnitude less than rank_atol are missing from S. For the cotangents, an arbitrary number of singular vectors or singular values can be missing, i.e. ΔU and ΔVᴴ can have sizes (m, pU) and (pV, n) respectively, whereas diagview(ΔS) can have length pS.

A warning will be printed if the cotangents are not gauge-invariant, i.e. if the anti-hermitian part of U' * ΔU + Vᴴ * ΔVᴴ', restricted to rows i and columns j for which abs(S[i] - S[j]) < degeneracy_atol, is not small compared to gauge_atol.

svd_full(A; kwargs...) -> U, S, Vᴴ
svd_full(A, alg::AbstractAlgorithm) -> U, S, Vᴴ
svd_full!(A, [USVᴴ]; kwargs...) -> U, S, Vᴴ
svd_full!(A, [USVᴴ], alg::AbstractAlgorithm) -> U, S, Vᴴ

Compute the full singular value decomposition (SVD) of the rectangular matrix A of size (m, n), such that A = U * S * Vᴴ. Here, U and Vᴴ are unitary matrices of size (m, m) and (n, n) respectively, and S is a diagonal matrix of size (m, n).

See also svd_compact(!), svd_vals(!) and svd_trunc(!).

svd_full(A; kwargs...) -> U, S, Vᴴ
svd_full(A, alg::AbstractAlgorithm) -> U, S, Vᴴ
svd_full!(A, [USVᴴ]; kwargs...) -> U, S, Vᴴ
svd_full!(A, [USVᴴ], alg::AbstractAlgorithm) -> U, S, Vᴴ

Compute the full singular value decomposition (SVD) of the rectangular matrix A of size (m, n), such that A = U * S * Vᴴ. Here, U and Vᴴ are unitary matrices of size (m, m) and (n, n) respectively, and S is a diagonal matrix of size (m, n).

See also svd_compact(!), svd_vals(!) and svd_trunc(!).

svd_trunc(A; kwargs...) -> U, S, Vᴴ
svd_trunc(A, alg::AbstractAlgorithm) -> U, S, Vᴴ
svd_trunc!(A, [USVᴴ]; kwargs...) -> U, S, Vᴴ
svd_trunc!(A, [USVᴴ], alg::AbstractAlgorithm) -> U, S, Vᴴ

Compute a partial or truncated singular value decomposition (SVD) of A, such that A * (Vᴴ)' = U * S. Here, U is an isometric matrix (orthonormal columns) of size (m, k), whereas Vᴴ is a matrix of size (k, n) with orthonormal rows and S is a square diagonal matrix of size (k, k), with k is set by the truncation strategy.

See also svd_full(!), svd_compact(!) and svd_vals(!).

svd_trunc(A; kwargs...) -> U, S, Vᴴ
svd_trunc(A, alg::AbstractAlgorithm) -> U, S, Vᴴ
svd_trunc!(A, [USVᴴ]; kwargs...) -> U, S, Vᴴ
svd_trunc!(A, [USVᴴ], alg::AbstractAlgorithm) -> U, S, Vᴴ

Compute a partial or truncated singular value decomposition (SVD) of A, such that A * (Vᴴ)' = U * S. Here, U is an isometric matrix (orthonormal columns) of size (m, k), whereas Vᴴ is a matrix of size (k, n) with orthonormal rows and S is a square diagonal matrix of size (k, k), with k is set by the truncation strategy.

See also svd_full(!), svd_compact(!) and svd_vals(!).

svd_vals(A; kwargs...) -> S
svd_vals(A, alg::AbstractAlgorithm) -> S
svd_vals!(A, [S]; kwargs...) -> S
svd_vals!(A, [S], alg::AbstractAlgorithm) -> S

Compute the vector of singular values of A, such that for an M×N matrix A, S is a vector of size K = min(M, N), the number of kept singular values.

See also svd_full(!), svd_compact(!) and svd_trunc(!).

svd_vals(A; kwargs...) -> S
svd_vals(A, alg::AbstractAlgorithm) -> S
svd_vals!(A, [S]; kwargs...) -> S
svd_vals!(A, [S], alg::AbstractAlgorithm) -> S

Compute the vector of singular values of A, such that for an M×N matrix A, S is a vector of size K = min(M, N), the number of kept singular values.

See also svd_full(!), svd_compact(!) and svd_trunc(!).

truncabove(atol::Real; by=abs)

Truncation strategy to discard the values that are larger than atol according to by.

truncate!(f, out, strategy::TruncationStrategy)

Generic interface for post-truncating a decomposition, specified in out.

truncrank(howmany::Int; by=abs, rev=true)

Truncation strategy to keep the first howmany values when sorted according to by or the last howmany if rev is true.

trunctol(atol::Real; by=abs)

Truncation strategy to discard the values that are smaller than atol according to by.

@algdef AlgorithmName

Convenience macro to define an algorithm AlgorithmName that accepts generic keywords. This defines an exported alias for Algorithm{:AlgorithmName} along with some utility methods.

@check_scalar(x, y, [op], [eltype])

Check if eltype(x) == op(eltype(y)) and throw an error if not. By default op = identity and `eltype = eltype'.

@check_size(x, sz, [size])

Check if size(x) == sz and throw an error if not. By default, size = size.

@functiondef f

Convenience macro to define the boilerplate code that dispatches between several versions of f and f!. By default, this enables the following signatures to be defined in terms of the final f!(A, out, alg::Algorithm).

    f(A; kwargs...)
    f(A, alg::Algorithm)
    f!(A, [out]; kwargs...)
    f!(A, alg::Algorithm)

See also copy_input, select_algorithm and initialize_output.