Sparse Linear Algebra (API)
Sparse Linear Algebra
Sparse matrix solvers call functions from SuiteSparse.
The following factorizations are available:
Type | Description |
---|---|
CHOLMOD.Factor | Cholesky and LDLt factorizations |
UMFPACK.UmfpackLU | LU factorization |
SPQR.QRSparse | QR factorization |
LinearAlgebra.cholesky
— Functioncholesky(A::SparseMatrixCSC; shift = 0.0, check = true, perm = nothing) -> CHOLMOD.Factor
Compute the Cholesky factorization of a sparse positive definite matrix A
. A
must be a SparseMatrixCSC
or a Symmetric
/Hermitian
view of a SparseMatrixCSC
. Note that even if A
doesn't have the type tag, it must still be symmetric or Hermitian. If perm
is not given, a fill-reducing permutation is used. F = cholesky(A)
is most frequently used to solve systems of equations with F\b
, but also the methods diag
, det
, and logdet
are defined for F
. You can also extract individual factors from F
, using F.L
. However, since pivoting is on by default, the factorization is internally represented as A == P'*L*L'*P
with a permutation matrix P
; using just L
without accounting for P
will give incorrect answers. To include the effects of permutation, it's typically preferable to extract "combined" factors like PtL = F.PtL
(the equivalent of P'*L
) and LtP = F.UP
(the equivalent of L'*P
).
When check = true
, an error is thrown if the decomposition fails. When check = false
, responsibility for checking the decomposition's validity (via issuccess
) lies with the user.
Setting the optional shift
keyword argument computes the factorization of A+shift*I
instead of A
. If the perm
argument is provided, it should be a permutation of 1:size(A,1)
giving the ordering to use (instead of CHOLMOD's default AMD ordering).
Examples
In the following example, the fill-reducing permutation used is [3, 2, 1]
. If perm
is set to 1:3
to enforce no permutation, the number of nonzero elements in the factor is 6.
julia> A = [2 1 1; 1 2 0; 1 0 2]
3×3 Matrix{Int64}:
2 1 1
1 2 0
1 0 2
julia> C = cholesky(sparse(A))
SparseArrays.CHOLMOD.Factor{Float64, Int64}
type: LLt
method: simplicial
maxnnz: 5
nnz: 5
success: true
julia> C.p
3-element Vector{Int64}:
3
2
1
julia> L = sparse(C.L);
julia> Matrix(L)
3×3 Matrix{Float64}:
1.41421 0.0 0.0
0.0 1.41421 0.0
0.707107 0.707107 1.0
julia> L * L' ≈ A[C.p, C.p]
true
julia> P = sparse(1:3, C.p, ones(3))
3×3 SparseMatrixCSC{Float64, Int64} with 3 stored entries:
⋅ ⋅ 1.0
⋅ 1.0 ⋅
1.0 ⋅ ⋅
julia> P' * L * L' * P ≈ A
true
julia> C = cholesky(sparse(A), perm=1:3)
SparseArrays.CHOLMOD.Factor{Float64, Int64}
type: LLt
method: simplicial
maxnnz: 6
nnz: 6
success: true
julia> L = sparse(C.L);
julia> Matrix(L)
3×3 Matrix{Float64}:
1.41421 0.0 0.0
0.707107 1.22474 0.0
0.707107 -0.408248 1.1547
julia> L * L' ≈ A
true
This method uses the CHOLMOD[ACM887][DavisHager2009] library from SuiteSparse. CHOLMOD only supports real or complex types in single or double precision. Input matrices not of those element types will be converted to these types as appropriate.
Many other functions from CHOLMOD are wrapped but not exported from the Base.SparseArrays.CHOLMOD
module.
LinearAlgebra.cholesky!
— Functioncholesky!(F::CHOLMOD.Factor, A::SparseMatrixCSC; shift = 0.0, check = true) -> CHOLMOD.Factor
Compute the Cholesky ($LL'$) factorization of A
, reusing the symbolic factorization F
. A
must be a SparseMatrixCSC
or a Symmetric
/ Hermitian
view of a SparseMatrixCSC
. Note that even if A
doesn't have the type tag, it must still be symmetric or Hermitian.
See also cholesky
.
This method uses the CHOLMOD library from SuiteSparse, which only supports real or complex types in single or double precision. Input matrices not of those element types will be converted to these types as appropriate.
LinearAlgebra.lowrankdowndate
— Functionlowrankdowndate(F::CHOLMOD.Factor, C::AbstractArray) -> FF::CHOLMOD.Factor
Get an LDLt
Factorization of A + C*C'
given an LDLt
or LLt
factorization F
of A
.
The returned factor is always an LDLt
factorization.
See also lowrankdowndate!
, lowrankupdate
, lowrankupdate!
.
LinearAlgebra.lowrankdowndate!
— Functionlowrankdowndate!(F::CHOLMOD.Factor, C::AbstractArray)
Update an LDLt
or LLt
Factorization F
of A
to a factorization of A - C*C'
.
LLt
factorizations are converted to LDLt
.
See also lowrankdowndate
, lowrankupdate
, lowrankupdate!
.
SparseArrays.CHOLMOD.lowrankupdowndate!
— Functionlowrankupdowndate!(F::CHOLMOD.Factor, C::Sparse, update::Cint)
Update an LDLt
or LLt
Factorization F
of A
to a factorization of A ± C*C'
.
If sparsity preserving factorization is used, i.e. L*L' == P*A*P'
then the new factor will be L*L' == P*A*P' + C'*C
update
: Cint(1)
for A + CC'
, Cint(0)
for A - CC'
LinearAlgebra.ldlt
— Functionldlt(A::SparseMatrixCSC; shift = 0.0, check = true, perm=nothing) -> CHOLMOD.Factor
Compute the $LDL'$ factorization of a sparse matrix A
. A
must be a SparseMatrixCSC
or a Symmetric
/Hermitian
view of a SparseMatrixCSC
. Note that even if A
doesn't have the type tag, it must still be symmetric or Hermitian. A fill-reducing permutation is used. F = ldlt(A)
is most frequently used to solve systems of equations A*x = b
with F\b
. The returned factorization object F
also supports the methods diag
, det
, logdet
, and inv
. You can extract individual factors from F
using F.L
. However, since pivoting is on by default, the factorization is internally represented as A == P'*L*D*L'*P
with a permutation matrix P
; using just L
without accounting for P
will give incorrect answers. To include the effects of permutation, it is typically preferable to extract "combined" factors like PtL = F.PtL
(the equivalent of P'*L
) and LtP = F.UP
(the equivalent of L'*P
). The complete list of supported factors is :L, :PtL, :D, :UP, :U, :LD, :DU, :PtLD, :DUP
.
When check = true
, an error is thrown if the decomposition fails. When check = false
, responsibility for checking the decomposition's validity (via issuccess
) lies with the user.
Setting the optional shift
keyword argument computes the factorization of A+shift*I
instead of A
. If the perm
argument is provided, it should be a permutation of 1:size(A,1)
giving the ordering to use (instead of CHOLMOD's default AMD ordering).
This method uses the CHOLMOD[ACM887][DavisHager2009] library from SuiteSparse. CHOLMOD only supports real or complex types in single or double precision. Input matrices not of those element types will be converted to these types as appropriate.
Many other functions from CHOLMOD are wrapped but not exported from the Base.SparseArrays.CHOLMOD
module.
LinearAlgebra.qr
— Functionqr(A::SparseMatrixCSC; tol=_default_tol(A), ordering=ORDERING_DEFAULT) -> QRSparse
Compute the QR
factorization of a sparse matrix A
. Fill-reducing row and column permutations are used such that F.R = F.Q'*A[F.prow,F.pcol]
. The main application of this type is to solve least squares or underdetermined problems with \
. The function calls the C library SPQR[ACM933].
qr(A::SparseMatrixCSC)
uses the SPQR library that is part of SuiteSparse. As this library only supports sparse matrices with Float64
or ComplexF64
elements, as of Julia v1.4 qr
converts A
into a copy that is of type SparseMatrixCSC{Float64}
or SparseMatrixCSC{ComplexF64}
as appropriate.
Examples
julia> A = sparse([1,2,3,4], [1,1,2,2], [1.0,1.0,1.0,1.0])
4×2 SparseMatrixCSC{Float64, Int64} with 4 stored entries:
1.0 ⋅
1.0 ⋅
⋅ 1.0
⋅ 1.0
julia> qr(A)
SparseArrays.SPQR.QRSparse{Float64, Int64}
Q factor:
4×4 SparseArrays.SPQR.QRSparseQ{Float64, Int64}
R factor:
2×2 SparseMatrixCSC{Float64, Int64} with 2 stored entries:
-1.41421 ⋅
⋅ -1.41421
Row permutation:
4-element Vector{Int64}:
1
3
4
2
Column permutation:
2-element Vector{Int64}:
1
2
LinearAlgebra.lu
— Functionlu(A::AbstractSparseMatrixCSC; check = true, q = nothing, control = get_umfpack_control()) -> F::UmfpackLU
Compute the LU factorization of a sparse matrix A
.
For sparse A
with real or complex element type, the return type of F
is UmfpackLU{Tv, Ti}
, with Tv
= Float64
or ComplexF64
respectively and Ti
is an integer type (Int32
or Int64
).
When check = true
, an error is thrown if the decomposition fails. When check = false
, responsibility for checking the decomposition's validity (via issuccess
) lies with the user.
The permutation q
can either be a permutation vector or nothing
. If no permutation vector is provided or q
is nothing
, UMFPACK's default is used. If the permutation is not zero-based, a zero-based copy is made.
The control
vector defaults to the Julia SparseArrays package's default configuration for UMFPACK (NB: this is modified from the UMFPACK defaults to disable iterative refinement), but can be changed by passing a vector of length UMFPACK_CONTROL
, see the UMFPACK manual for possible configurations. For example to reenable iterative refinement:
umfpack_control = SparseArrays.UMFPACK.get_umfpack_control(Float64, Int64) # read Julia default configuration for a Float64 sparse matrix
SparseArrays.UMFPACK.show_umf_ctrl(umfpack_control) # optional - display values
umfpack_control[SparseArrays.UMFPACK.JL_UMFPACK_IRSTEP] = 2.0 # reenable iterative refinement (2 is UMFPACK default max iterative refinement steps)
Alu = lu(A; control = umfpack_control)
x = Alu \ b # solve Ax = b, including UMFPACK iterative refinement
The individual components of the factorization F
can be accessed by indexing:
Component | Description |
---|---|
L | L (lower triangular) part of LU |
U | U (upper triangular) part of LU |
p | right permutation Vector |
q | left permutation Vector |
Rs | Vector of scaling factors |
: | (L,U,p,q,Rs) components |
The relation between F
and A
is
F.L*F.U == (F.Rs .* A)[F.p, F.q]
F
further supports the following functions:
See also lu!
lu(A::AbstractSparseMatrixCSC)
uses the UMFPACK[ACM832] library that is part of SuiteSparse. As this library only supports sparse matrices with Float64
or ComplexF64
elements, lu
converts A
into a copy that is of type SparseMatrixCSC{Float64}
or SparseMatrixCSC{ComplexF64}
as appropriate.
- ACM887Chen, Y., Davis, T. A., Hager, W. W., & Rajamanickam, S. (2008). Algorithm 887: CHOLMOD, Supernodal Sparse Cholesky Factorization and Update/Downdate. ACM Trans. Math. Softw., 35(3). doi:10.1145/1391989.1391995
- DavisHager2009Davis, Timothy A., & Hager, W. W. (2009). Dynamic Supernodes in Sparse Cholesky Update/Downdate and Triangular Solves. ACM Trans. Math. Softw., 35(4). doi:10.1145/1462173.1462176
- ACM933Foster, L. V., & Davis, T. A. (2013). Algorithm 933: Reliable Calculation of Numerical Rank, Null Space Bases, Pseudoinverse Solutions, and Basic Solutions Using SuitesparseQR. ACM Trans. Math. Softw., 40(1). doi:10.1145/2513109.2513116
- ACM832Davis, Timothy A. (2004b). Algorithm 832: UMFPACK V4.3–-an Unsymmetric-Pattern Multifrontal Method. ACM Trans. Math. Softw., 30(2), 196–199. doi:10.1145/992200.992206