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, sortby::Function, rev::Bool)

Truncation strategy to keep the first howmany values when sorted according to sortby or the last howmany if rev is 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.

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(!).

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.

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, atol::Real=0, rtol::Real=0, alg]) -> N
left_null!(A, [N]; [kind::Symbol, atol::Real=0, rtol::Real=0, alg]) -> 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 atol and rtol can be used to control the precision in determining the rank of A via its singular values.

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 == :qrpos: N is computed using the positive QR decomposition. This requires iszero(atol) && iszero(rtol) and left_null!(A, [N], kind=:qrpos) is equivalent to qr_null!(A, [N], alg) with a default value alg = select_algorithm(qr_compact!, A; positive=true)

• kind == :qr: N is computed using the (nonpositive) QR decomposition. This requires iszero(atol) && iszero(rtol) 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)

• 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 :qrpos when iszero(atol) && iszero(rtol) and :svd otherwise. Finally, finer control is obtained by providing an explicit algorithm using the alg keyword argument, which should be compatible with the chosen or default value of kind.

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

left_orth(A; [kind::Symbol, atol::Real=0, rtol::Real=0, alg]) -> V, C
left_orth!(A, [VC]; [kind::Symbol, atol::Real=0, rtol::Real=0, alg]) -> 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 atol and rtol can be used to control the precision in determining the rank of A via its singular values.

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

When kind is provided, its possible values are

• kind == :qrpos: V and C are computed using the positive QR decomposition. This requires iszero(atol) && iszero(rtol) 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 == :qr: V and C are computed using the QR decomposition, This requires iszero(atol) && iszero(rtol) and left_orth!(A, [VC]) is equivalent to qr_compact!(A, [VC], alg) with a default value alg = select_algorithm(qr_compact!, A)

• kind == :polar: V and C are computed using the polar decomposition, This requires iszero(atol) && iszero(rtol) 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!, where V will contain the left singular vectors corresponding to the singular values that are larger than max(atol, rtol * σ₁), where σ₁ is the largest singular value of A. C is computed as the product of the singular values and the right singular vectors, i.e. with U, S, Vᴴ = svd_trunc!(A), we have V = U and C = S * Vᴴ.

When kind is not provided, the default value is :qrpos when iszero(atol) && iszero(rtol) and :svd otherwise. Finally, finer control is obtained by providing an explicit algorithm using the alg keyword argument, which should be compatible with the chosen or default value of kind.

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, atol::Real=0, rtol::Real=0, alg]) -> Nᴴ
right_null!(A, [Nᴴ]; [kind::Symbol, atol::Real=0, rtol::Real=0, alg]) -> 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 atol and rtol can be used to control the precision in determining the rank of A via its singular values.

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 == :lqpos: Nᴴ is computed using the positive LQ decomposition. This requires iszero(atol) && iszero(rtol) and right_null!(A, [Nᴴ], kind=:lqpos) is equivalent to lq_null!(A, [Nᴴ], alg) with a default value alg = select_algorithm(lq_compact!, A; positive=true)

• kind == :lq: Nᴴ is computed using the (nonpositive) LQ decomposition. This requires iszero(atol) && iszero(rtol) 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)

• 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 :lqpos when iszero(atol) && iszero(rtol) and :svd otherwise. Finally, finer control is obtained by providing an explicit algorithm using the alg keyword argument, which should be compatible with the chosen or default value of kind.

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

right_orth(A; [kind::Symbol, atol::Real=0, rtol::Real=0, alg]) -> C, Vᴴ
right_orth!(A, [CVᴴ]; [kind::Symbol, atol::Real=0, rtol::Real=0, alg]) -> 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 atol and rtol can be used to control the precision in determining the rank of A via its singular values.

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

When kind is provided, its possible values are

• kind == :lqpos: C and Vᴴ are computed using the positive QR decomposition. This requires iszero(atol) && iszero(rtol) 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 == :lq: C and Vᴴ are computed using the QR decomposition, This requires iszero(atol) && iszero(rtol) and right_orth!(A, [CVᴴ]) is equivalent to lq_compact!(A, [CVᴴ], alg) with a default value alg = select_algorithm(lq_compact!, A))

• kind == :polar: C and Vᴴ are computed using the polar decomposition, This requires iszero(atol) && iszero(rtol) 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!, where V = adjoint(Vᴴ) will contain the right singular vectors corresponding to the singular values that are larger than max(atol, rtol * σ₁), where σ₁ is the largest singular value of A. C is computed as the product of the singular values and the right singular vectors, i.e. with U, S, Vᴴ = svd_trunc!(A), we have C = rmul!(U, S) and Vᴴ = Vᴴ.

When kind is not provided, the default value is :lqpos when iszero(atol) && iszero(rtol) and :svd otherwise. Finally, finer control is obtained by providing an explicit algorithm using the alg keyword argument, which should be compatible with the chosen or default value of kind.

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(!).

select_algorithm(f, A; kwargs...)

Given some keyword arguments and an input A, decide on an algrithm to use for implementing the function f on inputs of type A.

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)

Truncation strategy to discard the values that are larger than atol in absolute value.

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)

Truncation strategy to discard the values that are smaller than atol in absolute value.

@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.