Julia has support for sparse vectors andsparse matrices in theSparseArrays stdlib module. Sparse arrays are arrays that contain enough zeros that storing them in a special data structure leads to savings in space and execution time, compared to dense arrays.
External packages which implement different sparse storage types, multidimensional sparse arrays, and more can be found inNoteworthy External Sparse Packages
In Julia, sparse matrices are stored in theCompressed Sparse Column (CSC) format. Julia sparse matrices have the typeSparseMatrixCSC{Tv,Ti}, whereTv is the type of the stored values, andTi is the integer type for storing column pointers and row indices. The internal representation ofSparseMatrixCSC is as follows:
struct SparseMatrixCSC{Tv,Ti<:Integer} <: AbstractSparseMatrixCSC{Tv,Ti} m::Int # Number of rows n::Int # Number of columns colptr::Vector{Ti} # Column j is in colptr[j]:(colptr[j+1]-1) rowval::Vector{Ti} # Row indices of stored values nzval::Vector{Tv} # Stored values, typically nonzerosendThe compressed sparse column storage makes it easy and quick to access the elements in the column of a sparse matrix, whereas accessing the sparse matrix by rows is considerably slower. Operations such as insertion of previously unstored entries one at a time in the CSC structure tend to be slow. This is because all elements of the sparse matrix that are beyond the point of insertion have to be moved one place over.
All operations on sparse matrices are carefully implemented to exploit the CSC data structure for performance, and to avoid expensive operations.
If you have data in CSC format from a different application or library, and wish to import it in Julia, make sure that you use 1-based indexing. The row indices in every column need to be sorted, and if they are not, the matrix will display incorrectly. If yourSparseMatrixCSC object contains unsorted row indices, one quick way to sort them is by doing a double transpose. Since the transpose operation is lazy, make a copy to materialize each transpose.
In some applications, it is convenient to store explicit zero values in aSparseMatrixCSC. Theseare accepted by functions inBase (but there is no guarantee that they will be preserved in mutating operations). Such explicitly stored zeros are treated as structural nonzeros by many routines. Thennz function returns the number of elements explicitly stored in the sparse data structure, including non-structural zeros. In order to count the exact number of numerical nonzeros, usecount(!iszero, x), which inspects every stored element of a sparse matrix.dropzeros, and the in-placedropzeros!, can be used to remove stored zeros from the sparse matrix.
julia> A = sparse([1, 1, 2, 3], [1, 3, 2, 3], [0, 1, 2, 0])3×3 SparseMatrixCSC{Int64, Int64} with 4 stored entries: 0 ⋅ 1 ⋅ 2 ⋅ ⋅ ⋅ 0julia> dropzeros(A)3×3 SparseMatrixCSC{Int64, Int64} with 2 stored entries: ⋅ ⋅ 1 ⋅ 2 ⋅ ⋅ ⋅ ⋅Sparse vectors are stored in a close analog to compressed sparse column format for sparse matrices. In Julia, sparse vectors have the typeSparseVector{Tv,Ti} whereTv is the type of the stored values andTi the integer type for the indices. The internal representation is as follows:
struct SparseVector{Tv,Ti<:Integer} <: AbstractSparseVector{Tv,Ti} n::Int # Length of the sparse vector nzind::Vector{Ti} # Indices of stored values nzval::Vector{Tv} # Stored values, typically nonzerosendAs forSparseMatrixCSC, theSparseVector type can also contain explicitly stored zeros. (SeeSparse Matrix Storage.).
The simplest way to create a sparse array is to use a function equivalent to thezeros function that Julia provides for working with dense arrays. To produce a sparse array instead, you can use the same name with ansp prefix:
julia> spzeros(3)3-element SparseVector{Float64, Int64} with 0 stored entriesThesparse function is often a handy way to construct sparse arrays. For example, to construct a sparse matrix we can input a vectorI of row indices, a vectorJ of column indices, and a vectorV of stored values (this is also known as theCOO (coordinate) format).sparse(I,J,V) then constructs a sparse matrix such thatS[I[k], J[k]] = V[k]. The equivalent sparse vector constructor issparsevec, which takes the (row) index vectorI and the vectorV with the stored values and constructs a sparse vectorR such thatR[I[k]] = V[k].
julia> I = [1, 4, 3, 5]; J = [4, 7, 18, 9]; V = [1, 2, -5, 3];julia> S = sparse(I,J,V)5×18 SparseMatrixCSC{Int64, Int64} with 4 stored entries:⎡⠀⠈⠀⠀⠀⠀⠀⠀⢀⎤⎣⠀⠀⠀⠂⡀⠀⠀⠀⠀⎦julia> R = sparsevec(I,V)5-element SparseVector{Int64, Int64} with 4 stored entries: [1] = 1 [3] = -5 [4] = 2 [5] = 3The inverse of thesparse andsparsevec functions isfindnz, which retrieves the inputs used to create the sparse array (including stored entries equal to zero).findall(!iszero, x) returns the Cartesian indices of non-zero entries inx (not including stored entries equal to zero).
julia> findnz(S)([1, 4, 5, 3], [4, 7, 9, 18], [1, 2, 3, -5])julia> findall(!iszero, S)4-element Vector{CartesianIndex{2}}: CartesianIndex(1, 4) CartesianIndex(4, 7) CartesianIndex(5, 9) CartesianIndex(3, 18)julia> findnz(R)([1, 3, 4, 5], [1, -5, 2, 3])julia> findall(!iszero, R)4-element Vector{Int64}: 1 3 4 5Another way to create a sparse array is to convert a dense array into a sparse array using thesparse function:
julia> sparse(Matrix(1.0I, 5, 5))5×5 SparseMatrixCSC{Float64, Int64} with 5 stored entries: 1.0 ⋅ ⋅ ⋅ ⋅ ⋅ 1.0 ⋅ ⋅ ⋅ ⋅ ⋅ 1.0 ⋅ ⋅ ⋅ ⋅ ⋅ 1.0 ⋅ ⋅ ⋅ ⋅ ⋅ 1.0julia> sparse([1.0, 0.0, 1.0])3-element SparseVector{Float64, Int64} with 2 stored entries: [1] = 1.0 [3] = 1.0You can go in the other direction using theArray constructor. Theissparse function can be used to query if a matrix is sparse.
julia> issparse(spzeros(5))trueArithmetic operations on sparse matrices also work as they do on dense matrices. Indexing of, assignment into, and concatenation of sparse matrices work in the same way as dense matrices. Indexing operations, especially assignment, are expensive, when carried out one element at a time. In many cases it may be better to convert the sparse matrix into(I,J,V) format usingfindnz, manipulate the values or the structure in the dense vectors(I,J,V), and then reconstruct the sparse matrix.
The following table gives a correspondence between built-in methods on sparse matrices and their corresponding methods on dense matrix types. In general, methods that generate sparse matrices differ from their dense counterparts in that the resulting matrix follows the same sparsity pattern as a given sparse matrixS, or that the resulting sparse matrix has densityd, i.e. each matrix element has a probabilityd of being non-zero.
Details can be found in theSparse Vectors and Matrices section of the standard library reference.
| Sparse | Dense | Description |
|---|---|---|
spzeros(m,n) | zeros(m,n) | Creates am-by-n matrix of zeros. (spzeros(m,n) is empty.) |
sparse(I,n,n) | Matrix(I,n,n) | Creates an-by-n identity matrix. |
sparse(A) | Array(S) | Interconverts between dense and sparse formats. |
sprand(m,n,d) | rand(m,n) | Creates am-by-n random matrix (of densityd) with iid non-zero elements distributed uniformly on the half-open interval$[0, 1)$. |
sprandn(m,n,d) | randn(m,n) | Creates am-by-n random matrix (of densityd) with iid non-zero elements distributed according to the standard normal (Gaussian) distribution. |
sprandn(rng,m,n,d) | randn(rng,m,n) | Creates am-by-n random matrix (of densityd) with iid non-zero elements generated with therng random number generator |
Sparse matrix solvers call functions fromSuiteSparse. The following factorizations are available:
| Type | Description |
|---|---|
CHOLMOD.Factor | Cholesky and LDLt factorizations |
UMFPACK.UmfpackLU | LU factorization |
SPQR.QRSparse | QR factorization |
These factorizations are described in more detail in theSparse Linear Algebra API section:
SparseArrays.AbstractSparseArray —TypeAbstractSparseArray{Tv,Ti,N}Supertype forN-dimensional sparse arrays (or array-like types) with elements of typeTv and index typeTi.SparseMatrixCSC,SparseVector andSuiteSparse.CHOLMOD.Sparse are subtypes of this.
SparseArrays.AbstractSparseVector —TypeAbstractSparseVector{Tv,Ti}Supertype for one-dimensional sparse arrays (or array-like types) with elements of typeTv and index typeTi. Alias forAbstractSparseArray{Tv,Ti,1}.
SparseArrays.AbstractSparseMatrix —TypeAbstractSparseMatrix{Tv,Ti}Supertype for two-dimensional sparse arrays (or array-like types) with elements of typeTv and index typeTi. Alias forAbstractSparseArray{Tv,Ti,2}.
SparseArrays.SparseVector —TypeSparseVector{Tv,Ti<:Integer} <: AbstractSparseVector{Tv,Ti}Vector type for storing sparse vectors. Can be created by passing the length of the vector, asorted vector of non-zero indices, and a vector of non-zero values.
For instance, the vector[5, 6, 0, 7] can be represented as
SparseVector(4, [1, 2, 4], [5, 6, 7])This indicates that the element at index 1 is 5, at index 2 is 6, at index 3 iszero(Int), and at index 4 is 7.
It may be more convenient to create sparse vectors directly from dense vectors usingsparse as
sparse([5, 6, 0, 7])yields the same sparse vector.
SparseArrays.SparseMatrixCSC —TypeSparseMatrixCSC{Tv,Ti<:Integer} <: AbstractSparseMatrixCSC{Tv,Ti}Matrix type for storing sparse matrices in theCompressed Sparse Column format. The standard way of constructing SparseMatrixCSC is through thesparse function. See alsospzeros,spdiagm andsprand.
SparseArrays.sparse —Functionsparse(A)Convert an AbstractMatrixA into a sparse matrix.
Examples
julia> A = Matrix(1.0I, 3, 3)3×3 Matrix{Float64}: 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0julia> sparse(A)3×3 SparseMatrixCSC{Float64, Int64} with 3 stored entries: 1.0 ⋅ ⋅ ⋅ 1.0 ⋅ ⋅ ⋅ 1.0sparse(I, J, V,[ m, n, combine])Create a sparse matrixS of dimensionsm x n such thatS[I[k], J[k]] = V[k]. Thecombine function is used to combine duplicates. Ifm andn are not specified, they are set tomaximum(I) andmaximum(J) respectively. If thecombine function is not supplied,combine defaults to+ unless the elements ofV are Booleans in which casecombine defaults to|. All elements ofI must satisfy1 <= I[k] <= m, and all elements ofJ must satisfy1 <= J[k] <= n. Numerical zeros in (I,J,V) are retained as structural nonzeros; to drop numerical zeros, usedropzeros!.
For additional documentation and an expert driver, seeSparseArrays.sparse!.
Examples
julia> Is = [1; 2; 3];julia> Js = [1; 2; 3];julia> Vs = [1; 2; 3];julia> sparse(Is, Js, Vs)3×3 SparseMatrixCSC{Int64, Int64} with 3 stored entries: 1 ⋅ ⋅ ⋅ 2 ⋅ ⋅ ⋅ 3SparseArrays.sparse! —Functionsparse!(I::AbstractVector{Ti}, J::AbstractVector{Ti}, V::AbstractVector{Tv}, m::Integer, n::Integer, combine, klasttouch::Vector{Ti}, csrrowptr::Vector{Ti}, csrcolval::Vector{Ti}, csrnzval::Vector{Tv}, [csccolptr::Vector{Ti}], [cscrowval::Vector{Ti}, cscnzval::Vector{Tv}] ) where {Tv,Ti<:Integer}Parent of and expert driver forsparse; seesparse for basic usage. This method allows the user to provide preallocated storage forsparse's intermediate objects and result as described below. This capability enables more efficient successive construction ofSparseMatrixCSCs from coordinate representations, and also enables extraction of an unsorted-column representation of the result's transpose at no additional cost.
This method consists of three major steps: (1) Counting-sort the provided coordinate representation into an unsorted-row CSR form including repeated entries. (2) Sweep through the CSR form, simultaneously calculating the desired CSC form's column-pointer array, detecting repeated entries, and repacking the CSR form with repeated entries combined; this stage yields an unsorted-row CSR form with no repeated entries. (3) Counting-sort the preceding CSR form into a fully-sorted CSC form with no repeated entries.
Input arrayscsrrowptr,csrcolval, andcsrnzval constitute storage for the intermediate CSR forms and requirelength(csrrowptr) >= m + 1,length(csrcolval) >= length(I), andlength(csrnzval >= length(I)). Input arrayklasttouch, workspace for the second stage, requireslength(klasttouch) >= n. Optional input arrayscsccolptr,cscrowval, andcscnzval constitute storage for the returned CSC formS. If necessary, these are resized automatically to satisfylength(csccolptr) = n + 1,length(cscrowval) = nnz(S) andlength(cscnzval) = nnz(S); hence, ifnnz(S) is unknown at the outset, passing in empty vectors of the appropriate type (Vector{Ti}() andVector{Tv}() respectively) suffices, or calling thesparse! method neglectingcscrowval andcscnzval.
On return,csrrowptr,csrcolval, andcsrnzval contain an unsorted-column representation of the result's transpose.
You may reuse the input arrays' storage (I,J,V) for the output arrays (csccolptr,cscrowval,cscnzval). For example, you may callsparse!(I, J, V, csrrowptr, csrcolval, csrnzval, I, J, V). Note that they will be resized to satisfy the conditions above.
For the sake of efficiency, this method performs no argument checking beyond1 <= I[k] <= m and1 <= J[k] <= n. Use with care. Testing with--check-bounds=yes is wise.
This method runs inO(m, n, length(I)) time. The HALFPERM algorithm described in F. Gustavson, "Two fast algorithms for sparse matrices: multiplication and permuted transposition," ACM TOMS 4(3), 250-269 (1978) inspired this method's use of a pair of counting sorts.
SparseArrays.sparse!(I, J, V, [m, n, combine]) -> SparseMatrixCSCVariant ofsparse! that re-uses the input vectors (I,J,V) for the final matrix storage. After construction the input vectors will alias the matrix buffers;S.colptr === I,S.rowval === J, andS.nzval === V holds, and they will beresize!d as necessary.
Note that some work buffers will still be allocated. Specifically, this method is a convenience wrapper aroundsparse!(I, J, V, m, n, combine, klasttouch, csrrowptr, csrcolval, csrnzval, csccolptr, cscrowval, cscnzval) where this method allocatesklasttouch,csrrowptr,csrcolval, andcsrnzval of appropriate size, but reusesI,J, andV forcsccolptr,cscrowval, andcscnzval.
Argumentsm,n, andcombine defaults tomaximum(I),maximum(J), and+, respectively.
This method requires Julia version 1.10 or later.
SparseArrays.sparsevec —Functionsparsevec(I, V, [m, combine])Create a sparse vectorS of lengthm such thatS[I[k]] = V[k]. Duplicates are combined using thecombine function, which defaults to+ if nocombine argument is provided, unless the elements ofV are Booleans in which casecombine defaults to|.
Examples
julia> II = [1, 3, 3, 5]; V = [0.1, 0.2, 0.3, 0.2];julia> sparsevec(II, V)5-element SparseVector{Float64, Int64} with 3 stored entries: [1] = 0.1 [3] = 0.5 [5] = 0.2julia> sparsevec(II, V, 8, -)8-element SparseVector{Float64, Int64} with 3 stored entries: [1] = 0.1 [3] = -0.1 [5] = 0.2julia> sparsevec([1, 3, 1, 2, 2], [true, true, false, false, false])3-element SparseVector{Bool, Int64} with 3 stored entries: [1] = 1 [2] = 0 [3] = 1sparsevec(d::Dict, [m])Create a sparse vector of lengthm where the nonzero indices are keys from the dictionary, and the nonzero values are the values from the dictionary.
Examples
julia> sparsevec(Dict(1 => 3, 2 => 2))2-element SparseVector{Int64, Int64} with 2 stored entries: [1] = 3 [2] = 2sparsevec(A)Convert a vectorA into a sparse vector of lengthm.
Examples
julia> sparsevec([1.0, 2.0, 0.0, 0.0, 3.0, 0.0])6-element SparseVector{Float64, Int64} with 3 stored entries: [1] = 1.0 [2] = 2.0 [5] = 3.0Base.similar —Methodsimilar(A::AbstractSparseMatrixCSC{Tv,Ti}, [::Type{TvNew}, ::Type{TiNew}, m::Integer, n::Integer]) where {Tv,Ti}Create an uninitialized mutable array with the given element type, index type, and size, based upon the given sourceSparseMatrixCSC. The new sparse matrix maintains the structure of the original sparse matrix, except in the case where dimensions of the output matrix are different from the output.
The output matrix has zeros in the same locations as the input, but uninitialized values for the nonzero locations.
SparseArrays.issparse —Functionissparse(S)Returnstrue ifS is sparse, andfalse otherwise.
Examples
julia> sv = sparsevec([1, 4], [2.3, 2.2], 10)10-element SparseVector{Float64, Int64} with 2 stored entries: [1] = 2.3 [4] = 2.2julia> issparse(sv)truejulia> issparse(Array(sv))falseSparseArrays.nnz —Functionnnz(A)Returns the number of stored (filled) elements in a sparse array.
Examples
julia> A = sparse(2I, 3, 3)3×3 SparseMatrixCSC{Int64, Int64} with 3 stored entries: 2 ⋅ ⋅ ⋅ 2 ⋅ ⋅ ⋅ 2julia> nnz(A)3SparseArrays.findnz —Functionfindnz(A::SparseMatrixCSC)Return a tuple(I, J, V) whereI andJ are the row and column indices of the stored ("structurally non-zero") values in sparse matrixA, andV is a vector of the values.
Examples
julia> A = sparse([1 2 0; 0 0 3; 0 4 0])3×3 SparseMatrixCSC{Int64, Int64} with 4 stored entries: 1 2 ⋅ ⋅ ⋅ 3 ⋅ 4 ⋅julia> findnz(A)([1, 1, 3, 2], [1, 2, 2, 3], [1, 2, 4, 3])SparseArrays.spzeros —Functionspzeros([type,]m[,n])Create a sparse vector of lengthm or sparse matrix of sizem x n. This sparse array will not contain any nonzero values. No storage will be allocated for nonzero values during construction. The type defaults toFloat64 if not specified.
Examples
julia> spzeros(3, 3)3×3 SparseMatrixCSC{Float64, Int64} with 0 stored entries: ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅julia> spzeros(Float32, 4)4-element SparseVector{Float32, Int64} with 0 stored entriesspzeros([type], I::AbstractVector, J::AbstractVector, [m, n])Create a sparse matrixS of dimensionsm x n with structural zeros atS[I[k], J[k]].
This method can be used to construct the sparsity pattern of the matrix, and is more efficient than using e.g.sparse(I, J, zeros(length(I))).
For additional documentation and an expert driver, seeSparseArrays.spzeros!.
This methods requires Julia version 1.10 or later.
SparseArrays.spzeros! —Functionspzeros!(::Type{Tv}, I::AbstractVector{Ti}, J::AbstractVector{Ti}, m::Integer, n::Integer, klasttouch::Vector{Ti}, csrrowptr::Vector{Ti}, csrcolval::Vector{Ti}, [csccolptr::Vector{Ti}], [cscrowval::Vector{Ti}, cscnzval::Vector{Tv}]) where {Tv,Ti<:Integer}Parent of and expert driver forspzeros(I, J) allowing user to provide preallocated storage for intermediate objects. This method is tospzeros whatSparseArrays.sparse! is tosparse. See documentation forSparseArrays.sparse! for details and required buffer lengths.
This methods requires Julia version 1.10 or later.
SparseArrays.spzeros!(::Type{Tv}, I, J, [m, n]) -> SparseMatrixCSC{Tv}Variant ofspzeros! that re-uses the input vectorsI andJ for the final matrix storage. After construction the input vectors will alias the matrix buffers;S.colptr === I andS.rowval === J holds, and they will beresize!d as necessary.
Note that some work buffers will still be allocated. Specifically, this method is a convenience wrapper aroundspzeros!(Tv, I, J, m, n, klasttouch, csrrowptr, csrcolval, csccolptr, cscrowval) where this method allocatesklasttouch,csrrowptr, andcsrcolval of appropriate size, but reusesI andJ forcsccolptr andcscrowval.
Argumentsm andn defaults tomaximum(I) andmaximum(J).
This method requires Julia version 1.10 or later.
SparseArrays.spdiagm —Functionspdiagm(kv::Pair{<:Integer,<:AbstractVector}...)spdiagm(m::Integer, n::Integer, kv::Pair{<:Integer,<:AbstractVector}...)Construct a sparse diagonal matrix fromPairs of vectors and diagonals. Each vectorkv.second will be placed on thekv.first diagonal. By default, the matrix is square and its size is inferred fromkv, but a non-square sizem×n (padded with zeros as needed) can be specified by passingm,n as the first arguments.
Examples
julia> spdiagm(-1 => [1,2,3,4], 1 => [4,3,2,1])5×5 SparseMatrixCSC{Int64, Int64} with 8 stored entries: ⋅ 4 ⋅ ⋅ ⋅ 1 ⋅ 3 ⋅ ⋅ ⋅ 2 ⋅ 2 ⋅ ⋅ ⋅ 3 ⋅ 1 ⋅ ⋅ ⋅ 4 ⋅spdiagm(v::AbstractVector)spdiagm(m::Integer, n::Integer, v::AbstractVector)Construct a sparse matrix with elements of the vector as diagonal elements. By default (no givenm andn), the matrix is square and its size is given bylength(v), but a non-square sizem×n can be specified by passingm andn as the first arguments.
These functions require at least Julia 1.6.
Examples
julia> spdiagm([1,2,3])3×3 SparseMatrixCSC{Int64, Int64} with 3 stored entries: 1 ⋅ ⋅ ⋅ 2 ⋅ ⋅ ⋅ 3julia> spdiagm(sparse([1,0,3]))3×3 SparseMatrixCSC{Int64, Int64} with 2 stored entries: 1 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 3SparseArrays.sparse_hcat —Functionsparse_hcat(A...)Concatenate along dimension 2. Return a SparseMatrixCSC object.
This method was added in Julia 1.8. It mimics previous concatenation behavior, where the concatenation with specialized "sparse" matrix types from LinearAlgebra.jl automatically yielded sparse output even in the absence of any SparseArray argument.
SparseArrays.sparse_vcat —Functionsparse_vcat(A...)Concatenate along dimension 1. Return a SparseMatrixCSC object.
This method was added in Julia 1.8. It mimics previous concatenation behavior, where the concatenation with specialized "sparse" matrix types from LinearAlgebra.jl automatically yielded sparse output even in the absence of any SparseArray argument.
SparseArrays.sparse_hvcat —Functionsparse_hvcat(rows::Tuple{Vararg{Int}}, values...)Sparse horizontal and vertical concatenation in one call. This function is called for block matrix syntax. The first argument specifies the number of arguments to concatenate in each block row.
This method was added in Julia 1.8. It mimics previous concatenation behavior, where the concatenation with specialized "sparse" matrix types from LinearAlgebra.jl automatically yielded sparse output even in the absence of any SparseArray argument.
SparseArrays.blockdiag —Functionblockdiag(A...)Concatenate matrices block-diagonally. Currently only implemented for sparse matrices.
Examples
julia> blockdiag(sparse(2I, 3, 3), sparse(4I, 2, 2))5×5 SparseMatrixCSC{Int64, Int64} with 5 stored entries: 2 ⋅ ⋅ ⋅ ⋅ ⋅ 2 ⋅ ⋅ ⋅ ⋅ ⋅ 2 ⋅ ⋅ ⋅ ⋅ ⋅ 4 ⋅ ⋅ ⋅ ⋅ ⋅ 4SparseArrays.sprand —Functionsprand([rng],[T::Type],m,[n],p::AbstractFloat)sprand([rng],m,[n],p::AbstractFloat,[rfn=rand])Create a random lengthm sparse vector orm byn sparse matrix, in which the probability of any element being nonzero is independently given byp (and hence the mean density of nonzeros is also exactlyp). The optionalrng argument specifies a random number generator, seeRandom Numbers. The optionalT argument specifies the element type, which defaults toFloat64.
By default, nonzero values are sampled from a uniform distribution using therand function, i.e. byrand(T), orrand(rng, T) ifrng is supplied; for the defaultT=Float64, this corresponds to nonzero values sampled uniformly in[0,1).
You can sample nonzero values from a different distribution by passing a customrfn function instead ofrand. This should be a functionrfn(k) that returns an array ofk random numbers sampled from the desired distribution; alternatively, ifrng is supplied, it should instead be a functionrfn(rng, k).
Examples
julia> sprand(Bool, 2, 2, 0.5)2×2 SparseMatrixCSC{Bool, Int64} with 2 stored entries: 1 1 ⋅ ⋅julia> sprand(Float64, 3, 0.75)3-element SparseVector{Float64, Int64} with 2 stored entries: [1] = 0.795547 [2] = 0.49425SparseArrays.sprandn —Functionsprandn([rng][,Type],m[,n],p::AbstractFloat)Create a random sparse vector of lengthm or sparse matrix of sizem byn with the specified (independent) probabilityp of any entry being nonzero, where nonzero values are sampled from the normal distribution. The optionalrng argument specifies a random number generator, seeRandom Numbers.
Specifying the output element typeType requires at least Julia 1.1.
Examples
julia> sprandn(2, 2, 0.75)2×2 SparseMatrixCSC{Float64, Int64} with 3 stored entries: -1.20577 ⋅ 0.311817 -0.234641SparseArrays.nonzeros —Functionnonzeros(A)Return a vector of the structural nonzero values in sparse arrayA. This includes zeros that are explicitly stored in the sparse array. The returned vector points directly to the internal nonzero storage ofA, and any modifications to the returned vector will mutateA as well. Seerowvals andnzrange.
Examples
julia> A = sparse(2I, 3, 3)3×3 SparseMatrixCSC{Int64, Int64} with 3 stored entries: 2 ⋅ ⋅ ⋅ 2 ⋅ ⋅ ⋅ 2julia> nonzeros(A)3-element Vector{Int64}: 2 2 2SparseArrays.rowvals —Functionrowvals(A::AbstractSparseMatrixCSC)Return a vector of the row indices ofA. Any modifications to the returned vector will mutateA as well. Providing access to how the row indices are stored internally can be useful in conjunction with iterating over structural nonzero values. See alsononzeros andnzrange.
Examples
julia> A = sparse(2I, 3, 3)3×3 SparseMatrixCSC{Int64, Int64} with 3 stored entries: 2 ⋅ ⋅ ⋅ 2 ⋅ ⋅ ⋅ 2julia> rowvals(A)3-element Vector{Int64}: 1 2 3SparseArrays.nzrange —Functionnzrange(A::AbstractSparseMatrixCSC, col::Integer)Return the range of indices to the structural nonzero values of a sparse matrix column. In conjunction withnonzeros androwvals, this allows for convenient iterating over a sparse matrix :
A = sparse(I,J,V)rows = rowvals(A)vals = nonzeros(A)m, n = size(A)for j = 1:n for i in nzrange(A, j) row = rows[i] val = vals[i] # perform sparse wizardry... endendAdding or removing nonzero elements to the matrix may invalidate thenzrange, one should not mutate the matrix while iterating.
nzrange(x::SparseVectorUnion, col)Give the range of indices to the structural nonzero values of a sparse vector. The column indexcol is ignored (assumed to be1).
SparseArrays.droptol! —Functiondroptol!(A::AbstractSparseMatrixCSC, tol)Removes stored values fromA whose absolute value is less than or equal totol.
droptol!(x::AbstractCompressedVector, tol)Removes stored values fromx whose absolute value is less than or equal totol.
SparseArrays.dropzeros! —Functiondropzeros!(x::AbstractCompressedVector)Removes stored numerical zeros fromx.
For an out-of-place version, seedropzeros. For algorithmic information, seefkeep!.
SparseArrays.dropzeros —Functiondropzeros(A::AbstractSparseMatrixCSC;)Generates a copy ofA and removes stored numerical zeros from that copy.
For an in-place version and algorithmic information, seedropzeros!.
Examples
julia> A = sparse([1, 2, 3], [1, 2, 3], [1.0, 0.0, 1.0])3×3 SparseMatrixCSC{Float64, Int64} with 3 stored entries: 1.0 ⋅ ⋅ ⋅ 0.0 ⋅ ⋅ ⋅ 1.0julia> dropzeros(A)3×3 SparseMatrixCSC{Float64, Int64} with 2 stored entries: 1.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1.0dropzeros(x::AbstractCompressedVector)Generates a copy ofx and removes numerical zeros from that copy.
For an in-place version and algorithmic information, seedropzeros!.
Examples
julia> A = sparsevec([1, 2, 3], [1.0, 0.0, 1.0])3-element SparseVector{Float64, Int64} with 3 stored entries: [1] = 1.0 [2] = 0.0 [3] = 1.0julia> dropzeros(A)3-element SparseVector{Float64, Int64} with 2 stored entries: [1] = 1.0 [3] = 1.0SparseArrays.permute —Functionpermute(A::AbstractSparseMatrixCSC{Tv,Ti}, p::AbstractVector{<:Integer}, q::AbstractVector{<:Integer}) where {Tv,Ti}Bilaterally permuteA, returningPAQ (A[p,q]). Column-permutationq's length must matchA's column count (length(q) == size(A, 2)). Row-permutationp's length must matchA's row count (length(p) == size(A, 1)).
For expert drivers and additional information, seepermute!.
Examples
julia> A = spdiagm(0 => [1, 2, 3, 4], 1 => [5, 6, 7])4×4 SparseMatrixCSC{Int64, Int64} with 7 stored entries: 1 5 ⋅ ⋅ ⋅ 2 6 ⋅ ⋅ ⋅ 3 7 ⋅ ⋅ ⋅ 4julia> permute(A, [4, 3, 2, 1], [1, 2, 3, 4])4×4 SparseMatrixCSC{Int64, Int64} with 7 stored entries: ⋅ ⋅ ⋅ 4 ⋅ ⋅ 3 7 ⋅ 2 6 ⋅ 1 5 ⋅ ⋅julia> permute(A, [1, 2, 3, 4], [4, 3, 2, 1])4×4 SparseMatrixCSC{Int64, Int64} with 7 stored entries: ⋅ ⋅ 5 1 ⋅ 6 2 ⋅ 7 3 ⋅ ⋅ 4 ⋅ ⋅ ⋅Base.permute! —Methodpermute!(X::AbstractSparseMatrixCSC{Tv,Ti}, A::AbstractSparseMatrixCSC{Tv,Ti}, p::AbstractVector{<:Integer}, q::AbstractVector{<:Integer}, [C::AbstractSparseMatrixCSC{Tv,Ti}]) where {Tv,Ti}Bilaterally permuteA, storing resultPAQ (A[p,q]) inX. Stores intermediate result(AQ)^T (transpose(A[:,q])) in optional argumentC if present. Requires that none ofX,A, and, if present,C alias each other; to store resultPAQ back intoA, use the following method lackingX:
permute!(A::AbstractSparseMatrixCSC{Tv,Ti}, p::AbstractVector{<:Integer}, q::AbstractVector{<:Integer}[, C::AbstractSparseMatrixCSC{Tv,Ti}, [workcolptr::Vector{Ti}]]) where {Tv,Ti}X's dimensions must match those ofA (size(X, 1) == size(A, 1) andsize(X, 2) == size(A, 2)), andX must have enough storage to accommodate all allocated entries inA (length(rowvals(X)) >= nnz(A) andlength(nonzeros(X)) >= nnz(A)). Column-permutationq's length must matchA's column count (length(q) == size(A, 2)). Row-permutationp's length must matchA's row count (length(p) == size(A, 1)).
C's dimensions must match those oftranspose(A) (size(C, 1) == size(A, 2) andsize(C, 2) == size(A, 1)), andC must have enough storage to accommodate all allocated entries inA (length(rowvals(C)) >= nnz(A) andlength(nonzeros(C)) >= nnz(A)).
For additional (algorithmic) information, and for versions of these methods that forgo argument checking, see (unexported) parent methodsunchecked_noalias_permute! andunchecked_aliasing_permute!.
See alsopermute.
SparseArrays.halfperm! —Functionhalfperm!(X::AbstractSparseMatrixCSC{Tv,Ti}, A::AbstractSparseMatrixCSC{TvA,Ti}, q::AbstractVector{<:Integer}, f::Function = identity) where {Tv,TvA,Ti}Column-permute and transposeA, simultaneously applyingf to each entry ofA, storing the result(f(A)Q)^T (map(f, transpose(A[:,q]))) inX.
Element typeTv ofX must matchf(::TvA), whereTvA is the element type ofA.X's dimensions must match those oftranspose(A) (size(X, 1) == size(A, 2) andsize(X, 2) == size(A, 1)), andX must have enough storage to accommodate all allocated entries inA (length(rowvals(X)) >= nnz(A) andlength(nonzeros(X)) >= nnz(A)). Column-permutationq's length must matchA's column count (length(q) == size(A, 2)).
This method is the parent of several methods performing transposition and permutation operations onSparseMatrixCSCs. As this method performs no argument checking, prefer the safer child methods ([c]transpose[!],permute[!]) to direct use.
This method implements theHALFPERM algorithm described in F. Gustavson, "Two fast algorithms for sparse matrices: multiplication and permuted transposition," ACM TOMS 4(3), 250-269 (1978). The algorithm runs inO(size(A, 1), size(A, 2), nnz(A)) time and requires no space beyond that passed in.
SparseArrays.ftranspose! —Functionftranspose!(X::AbstractSparseMatrixCSC{Tv,Ti}, A::AbstractSparseMatrixCSC{Tv,Ti}, f::Function) where {Tv,Ti}TransposeA and store it inX while applying the functionf to the non-zero elements. Does not remove the zeros created byf.size(X) must be equal tosize(transpose(A)). No additional memory is allocated other than resizing the rowval and nzval ofX, if needed.
Seehalfperm!
LinearAlgebra.cholesky —Functioncholesky(A, NoPivot(); check = true) -> CholeskyCompute the Cholesky factorization of a dense symmetric positive definite matrixA and return aCholesky factorization. The matrixA can either be aSymmetric orHermitianAbstractMatrix or aperfectly symmetric or HermitianAbstractMatrix.
The triangular Cholesky factor can be obtained from the factorizationF viaF.L andF.U, whereA ≈ F.U' * F.U ≈ F.L * F.L'.
The following functions are available forCholesky objects:size,\,inv,det,logdet andisposdef.
If you have a matrixA that is slightly non-Hermitian due to roundoff errors in its construction, wrap it inHermitian(A) before passing it tocholesky in order to treat it as perfectly Hermitian.
Whencheck = true, an error is thrown if the decomposition fails. Whencheck = false, responsibility for checking the decomposition's validity (viaissuccess) lies with the user.
Examples
julia> A = [4. 12. -16.; 12. 37. -43.; -16. -43. 98.]3×3 Matrix{Float64}: 4.0 12.0 -16.0 12.0 37.0 -43.0 -16.0 -43.0 98.0julia> C = cholesky(A)Cholesky{Float64, Matrix{Float64}}U factor:3×3 UpperTriangular{Float64, Matrix{Float64}}: 2.0 6.0 -8.0 ⋅ 1.0 5.0 ⋅ ⋅ 3.0julia> C.U3×3 UpperTriangular{Float64, Matrix{Float64}}: 2.0 6.0 -8.0 ⋅ 1.0 5.0 ⋅ ⋅ 3.0julia> C.L3×3 LowerTriangular{Float64, Matrix{Float64}}: 2.0 ⋅ ⋅ 6.0 1.0 ⋅ -8.0 5.0 3.0julia> C.L * C.U == Atruecholesky(A, RowMaximum(); tol = 0.0, check = true) -> CholeskyPivotedCompute the pivoted Cholesky factorization of a dense symmetric positive semi-definite matrixA and return aCholeskyPivoted factorization. The matrixA can either be aSymmetric orHermitianAbstractMatrix or aperfectly symmetric or HermitianAbstractMatrix.
The triangular Cholesky factor can be obtained from the factorizationF viaF.L andF.U, and the permutation viaF.p, whereA[F.p, F.p] ≈ Ur' * Ur ≈ Lr * Lr' withUr = F.U[1:F.rank, :] andLr = F.L[:, 1:F.rank], or alternativelyA ≈ Up' * Up ≈ Lp * Lp' withUp = F.U[1:F.rank, invperm(F.p)] andLp = F.L[invperm(F.p), 1:F.rank].
The following functions are available forCholeskyPivoted objects:size,\,inv,det, andrank.
The argumenttol determines the tolerance for determining the rank. For negative values, the tolerance is the machine precision.
If you have a matrixA that is slightly non-Hermitian due to roundoff errors in its construction, wrap it inHermitian(A) before passing it tocholesky in order to treat it as perfectly Hermitian.
Whencheck = true, an error is thrown if the decomposition fails. Whencheck = false, responsibility for checking the decomposition's validity (viaissuccess) lies with the user.
Examples
julia> X = [1.0, 2.0, 3.0, 4.0];julia> A = X * X';julia> C = cholesky(A, RowMaximum(), check = false)CholeskyPivoted{Float64, Matrix{Float64}, Vector{Int64}}U factor with rank 1:4×4 UpperTriangular{Float64, Matrix{Float64}}: 4.0 2.0 3.0 1.0 ⋅ 0.0 6.0 2.0 ⋅ ⋅ 9.0 3.0 ⋅ ⋅ ⋅ 1.0permutation:4-element Vector{Int64}: 4 2 3 1julia> C.U[1:C.rank, :]' * C.U[1:C.rank, :] ≈ A[C.p, C.p]truejulia> l, u = C; # destructuring via iterationjulia> l == C.L && u == C.Utruecholesky(A::SparseMatrixCSC; shift = 0.0, check = true, perm = nothing) -> CHOLMOD.FactorCompute the Cholesky factorization of a sparse positive definite matrixA.A must be aSparseMatrixCSC or aSymmetric/Hermitian view of aSparseMatrixCSC. Note that even ifA doesn't have the type tag, it must still be symmetric or Hermitian. Ifperm is not given, a fill-reducing permutation is used.F = cholesky(A) is most frequently used to solve systems of equations withF\b, but also the methodsdiag,det, andlogdet are defined forF. You can also extract individual factors fromF, usingF.L. However, since pivoting is on by default, the factorization is internally represented asA == P'*L*L'*P with a permutation matrixP; using justL without accounting forP will give incorrect answers. To include the effects of permutation, it's typically preferable to extract "combined" factors likePtL = F.PtL (the equivalent ofP'*L) andLtP = F.UP (the equivalent ofL'*P).
Whencheck = true, an error is thrown if the decomposition fails. Whencheck = false, responsibility for checking the decomposition's validity (viaissuccess) lies with the user.
Setting the optionalshift keyword argument computes the factorization ofA+shift*I instead ofA. If theperm argument is provided, it should be a permutation of1: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]. Ifperm is set to1: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 2julia> C = cholesky(sparse(A))SparseArrays.CHOLMOD.Factor{Float64, Int64}type: LLtmethod: simplicialmaxnnz: 5nnz: 5success: truejulia> C.p3-element Vector{Int64}: 3 2 1julia> 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.0julia> L * L' ≈ A[C.p, C.p]truejulia> 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 ≈ Atruejulia> C = cholesky(sparse(A), perm=1:3)SparseArrays.CHOLMOD.Factor{Float64, Int64}type: LLtmethod: simplicialmaxnnz: 6nnz: 6success: truejulia> 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.1547julia> L * L' ≈ AtrueThis method uses the CHOLMOD[ACM887][DavisHager2009] library fromSuiteSparse. 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 theBase.SparseArrays.CHOLMOD module.
LinearAlgebra.cholesky! —Functioncholesky!(A::AbstractMatrix, NoPivot(); check = true) -> CholeskyThe same ascholesky, but saves space by overwriting the inputA, instead of creating a copy. AnInexactError exception is thrown if the factorization produces a number not representable by the element type ofA, e.g. for integer types.
Examples
julia> A = [1 2; 2 50]2×2 Matrix{Int64}: 1 2 2 50julia> cholesky!(A)ERROR: InexactError: Int64(6.782329983125268)Stacktrace:[...]cholesky!(A::AbstractMatrix, RowMaximum(); tol = 0.0, check = true) -> CholeskyPivotedThe same ascholesky, but saves space by overwriting the inputA, instead of creating a copy. AnInexactError exception is thrown if the factorization produces a number not representable by the element type ofA, e.g. for integer types.
cholesky!(F::CHOLMOD.Factor, A::SparseMatrixCSC; shift = 0.0, check = true) -> CHOLMOD.FactorCompute the Cholesky ($LL'$) factorization ofA, reusing the symbolic factorizationF.A must be aSparseMatrixCSC or aSymmetric/Hermitian view of aSparseMatrixCSC. Note that even ifA doesn't have the type tag, it must still be symmetric or Hermitian.
See alsocholesky.
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.lowrankupdate —Functionlowrankupdate(C::Cholesky, v::AbstractVector) -> CC::CholeskyUpdate a Cholesky factorizationC with the vectorv. IfA = C.U'C.U thenCC = cholesky(C.U'C.U + v*v') but the computation ofCC only usesO(n^2) operations.
lowrankupdate(F::CHOLMOD.Factor, C::AbstractArray) -> FF::CHOLMOD.FactorGet anLDLt Factorization ofA + C*C' given anLDLt orLLt factorizationF ofA.
The returned factor is always anLDLt factorization.
LinearAlgebra.lowrankupdate! —Functionlowrankupdate!(C::Cholesky, v::AbstractVector) -> CC::CholeskyUpdate a Cholesky factorizationC with the vectorv. IfA = C.U'C.U thenCC = cholesky(C.U'C.U + v*v') but the computation ofCC only usesO(n^2) operations. The input factorizationC is updated in place such that on exitC == CC. The vectorv is destroyed during the computation.
lowrankupdate!(F::CHOLMOD.Factor, C::AbstractArray)Update anLDLt orLLt FactorizationF ofA to a factorization ofA + C*C'.
LLt factorizations are converted toLDLt.
See alsolowrankupdate,lowrankdowndate,lowrankdowndate!.
LinearAlgebra.lowrankdowndate —Functionlowrankdowndate(C::Cholesky, v::AbstractVector) -> CC::CholeskyDowndate a Cholesky factorizationC with the vectorv. IfA = C.U'C.U thenCC = cholesky(C.U'C.U - v*v') but the computation ofCC only usesO(n^2) operations.
lowrankdowndate(F::CHOLMOD.Factor, C::AbstractArray) -> FF::CHOLMOD.FactorGet anLDLt Factorization ofA + C*C' given anLDLt orLLt factorizationF ofA.
The returned factor is always anLDLt factorization.
See alsolowrankdowndate!,lowrankupdate,lowrankupdate!.
LinearAlgebra.lowrankdowndate! —Functionlowrankdowndate!(C::Cholesky, v::AbstractVector) -> CC::CholeskyDowndate a Cholesky factorizationC with the vectorv. IfA = C.U'C.U thenCC = cholesky(C.U'C.U - v*v') but the computation ofCC only usesO(n^2) operations. The input factorizationC is updated in place such that on exitC == CC. The vectorv is destroyed during the computation.
lowrankdowndate!(F::CHOLMOD.Factor, C::AbstractArray)Update anLDLt orLLt FactorizationF ofA to a factorization ofA - C*C'.
LLt factorizations are converted toLDLt.
See alsolowrankdowndate,lowrankupdate,lowrankupdate!.
SparseArrays.CHOLMOD.lowrankupdowndate! —Functionlowrankupdowndate!(F::CHOLMOD.Factor, C::Sparse, update::Cint)Update anLDLt orLLt FactorizationF ofA to a factorization ofA ± C*C'.
If sparsity preserving factorization is used, i.e.L*L' == P*A*P' then the new factor will beL*L' == P*A*P' + C'*C
update:Cint(1) forA + CC',Cint(0) forA - CC'
LinearAlgebra.ldlt —Functionldlt(S::SymTridiagonal) -> LDLtCompute anLDLt (i.e.,$LDL^T$) factorization of the real symmetric tridiagonal matrixS such thatS = L*Diagonal(d)*L' whereL is a unit lower triangular matrix andd is a vector. The main use of anLDLt factorizationF = ldlt(S) is to solve the linear system of equationsSx = b withF\b.
See alsobunchkaufman for a similar, but pivoted, factorization of arbitrary symmetric or Hermitian matrices.
Examples
julia> S = SymTridiagonal([3., 4., 5.], [1., 2.])3×3 SymTridiagonal{Float64, Vector{Float64}}: 3.0 1.0 ⋅ 1.0 4.0 2.0 ⋅ 2.0 5.0julia> ldltS = ldlt(S);julia> b = [6., 7., 8.];julia> ldltS \ b3-element Vector{Float64}: 1.7906976744186047 0.627906976744186 1.3488372093023255julia> S \ b3-element Vector{Float64}: 1.7906976744186047 0.627906976744186 1.3488372093023255ldlt(A::SparseMatrixCSC; shift = 0.0, check = true, perm=nothing) -> CHOLMOD.FactorCompute the$LDL'$ factorization of a sparse matrixA.A must be aSparseMatrixCSC or aSymmetric/Hermitian view of aSparseMatrixCSC. Note that even ifA 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 equationsA*x = b withF\b. The returned factorization objectF also supports the methodsdiag,det,logdet, andinv. You can extract individual factors fromF usingF.L. However, since pivoting is on by default, the factorization is internally represented asA == P'*L*D*L'*P with a permutation matrixP; using justL without accounting forP will give incorrect answers. To include the effects of permutation, it is typically preferable to extract "combined" factors likePtL = F.PtL (the equivalent ofP'*L) andLtP = F.UP (the equivalent ofL'*P). The complete list of supported factors is:L, :PtL, :D, :UP, :U, :LD, :DU, :PtLD, :DUP.
Whencheck = true, an error is thrown if the decomposition fails. Whencheck = false, responsibility for checking the decomposition's validity (viaissuccess) lies with the user.
Setting the optionalshift keyword argument computes the factorization ofA+shift*I instead ofA. If theperm argument is provided, it should be a permutation of1:size(A,1) giving the ordering to use (instead of CHOLMOD's default AMD ordering).
This method uses the CHOLMOD[ACM887][DavisHager2009] library fromSuiteSparse. 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 theBase.SparseArrays.CHOLMOD module.
LinearAlgebra.lu —Functionlu(A::AbstractSparseMatrixCSC; check = true, q = nothing, control = get_umfpack_control()) -> F::UmfpackLUCompute the LU factorization of a sparse matrixA.
For sparseA with real or complex element type, the return type ofF isUmfpackLU{Tv, Ti}, withTv =Float64 orComplexF64 respectively andTi is an integer type (Int32 orInt64).
Whencheck = true, an error is thrown if the decomposition fails. Whencheck = false, responsibility for checking the decomposition's validity (viaissuccess) lies with the user.
The permutationq can either be a permutation vector ornothing. If no permutation vector is provided orq isnothing, UMFPACK's default is used. If the permutation is not zero-based, a zero-based copy is made.
Thecontrol 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 lengthUMFPACK_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 matrixSparseArrays.UMFPACK.show_umf_ctrl(umfpack_control) # optional - display valuesumfpack_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 refinementThe individual components of the factorizationF can be accessed by indexing:
| Component | Description |
|---|---|
L | L (lower triangular) part ofLU |
U | U (upper triangular) part ofLU |
p | right permutationVector |
q | left permutationVector |
Rs | Vector of scaling factors |
: | (L,U,p,q,Rs) components |
The relation betweenF andA is
F.L*F.U == (F.Rs .* A)[F.p, F.q]
F further supports the following functions:
See alsolu!
lu(A::AbstractSparseMatrixCSC) uses the UMFPACK[ACM832] library that is part ofSuiteSparse. As this library only supports sparse matrices withFloat64 orComplexF64 elements,lu convertsA into a copy that is of typeSparseMatrixCSC{Float64} orSparseMatrixCSC{ComplexF64} as appropriate.
lu(A, pivot = RowMaximum(); check = true, allowsingular = false) -> F::LUCompute the LU factorization ofA.
Whencheck = true, an error is thrown if the decomposition fails. Whencheck = false, responsibility for checking the decomposition's validity (viaissuccess) lies with the user.
By default, withcheck = true, an error is also thrown when the decomposition produces valid factors, but the upper-triangular factorU is rank-deficient. This may be changed by passingallowsingular = true.
In most cases, ifA is a subtypeS ofAbstractMatrix{T} with an element typeT supporting+,-,* and/, the return type isLU{T,S{T}}.
In general, LU factorization involves a permutation of the rows of the matrix (corresponding to theF.p output described below), known as "pivoting" (because it corresponds to choosing which row contains the "pivot", the diagonal entry ofF.U). One of the following pivoting strategies can be selected via the optionalpivot argument:
RowMaximum() (default): the standard pivoting strategy; the pivot corresponds to the element of maximum absolute value among the remaining, to be factorized rows. This pivoting strategy requires the element type to also supportabs and<. (This is generally the only numerically stable option for floating-point matrices.)RowNonZero(): the pivot corresponds to the first non-zero element among the remaining, to be factorized rows. (This corresponds to the typical choice in hand calculations, and is also useful for more general algebraic number types that supportiszero but notabs or<.)NoPivot(): pivoting turned off (will fail if a zero entry is encountered in a pivot position, even whenallowsingular = true).The individual components of the factorizationF can be accessed viagetproperty:
| Component | Description |
|---|---|
F.L | L (lower triangular) part ofLU |
F.U | U (upper triangular) part ofLU |
F.p | (right) permutationVector |
F.P | (right) permutationMatrix |
Iterating the factorization produces the componentsF.L,F.U, andF.p.
The relationship betweenF andA is
F.L*F.U == A[F.p, :]
F further supports the following functions:
| Supported function | LU | LU{T,Tridiagonal{T}} |
|---|---|---|
/ | ✓ | |
\ | ✓ | ✓ |
inv | ✓ | ✓ |
det | ✓ | ✓ |
logdet | ✓ | ✓ |
logabsdet | ✓ | ✓ |
size | ✓ | ✓ |
Theallowsingular keyword argument was added in Julia 1.11.
Examples
julia> A = [4 3; 6 3]2×2 Matrix{Int64}: 4 3 6 3julia> F = lu(A)LU{Float64, Matrix{Float64}, Vector{Int64}}L factor:2×2 Matrix{Float64}: 1.0 0.0 0.666667 1.0U factor:2×2 Matrix{Float64}: 6.0 3.0 0.0 1.0julia> F.L * F.U == A[F.p, :]truejulia> l, u, p = lu(A); # destructuring via iterationjulia> l == F.L && u == F.U && p == F.ptruejulia> lu([1 2; 1 2], allowsingular = true)LU{Float64, Matrix{Float64}, Vector{Int64}}L factor:2×2 Matrix{Float64}: 1.0 0.0 1.0 1.0U factor (rank-deficient):2×2 Matrix{Float64}: 1.0 2.0 0.0 0.0LinearAlgebra.qr —Functionqr(A::SparseMatrixCSC; tol=_default_tol(A), ordering=ORDERING_DEFAULT) -> QRSparseCompute theQR factorization of a sparse matrixA. Fill-reducing row and column permutations are used such thatF.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 ofSuiteSparse. As this library only supports sparse matrices withFloat64 orComplexF64 elements, as of Julia v1.4qr convertsA into a copy that is of typeSparseMatrixCSC{Float64} orSparseMatrixCSC{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.0julia> 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.41421Row permutation:4-element Vector{Int64}: 1 3 4 2Column permutation:2-element Vector{Int64}: 1 2qr(A, pivot = NoPivot(); blocksize) -> FCompute the QR factorization of the matrixA: an orthogonal (or unitary ifA is complex-valued) matrixQ, and an upper triangular matrixR such that
\[A = Q R\]
The returned objectF stores the factorization in a packed format:
ifpivot == ColumnNorm() thenF is aQRPivoted object,
otherwise if the element type ofA is a BLAS type (Float32,Float64,ComplexF32 orComplexF64), thenF is aQRCompactWY object,
otherwiseF is aQR object.
The individual components of the decompositionF can be retrieved via property accessors:
F.Q: the orthogonal/unitary matrixQF.R: the upper triangular matrixRF.p: the permutation vector of the pivot (QRPivoted only)F.P: the permutation matrix of the pivot (QRPivoted only)Each reference to the upper triangular factor viaF.R allocates a new array. It is therefore advisable to cache that array, say, byR = F.R and continue working withR.
Iterating the decomposition produces the componentsQ,R, and if extantp.
The following functions are available for theQR objects:inv,size, and\. WhenA is rectangular,\ will return a least squares solution and if the solution is not unique, the one with smallest norm is returned. WhenA is not full rank, factorization with (column) pivoting is required to obtain a minimum norm solution.
Multiplication with respect to either full/square or non-full/squareQ is allowed, i.e. bothF.Q*F.R andF.Q*A are supported. AQ matrix can be converted into a regular matrix withMatrix. This operation returns the "thin" Q factor, i.e., ifA ism×n withm>=n, thenMatrix(F.Q) yields anm×n matrix with orthonormal columns. To retrieve the "full" Q factor, anm×m orthogonal matrix, useF.Q*I orcollect(F.Q). Ifm<=n, thenMatrix(F.Q) yields anm×m orthogonal matrix.
The block size for QR decomposition can be specified by keyword argumentblocksize :: Integer whenpivot == NoPivot() andA isa StridedMatrix{<:BlasFloat}. It is ignored whenblocksize > minimum(size(A)). SeeQRCompactWY.
Theblocksize keyword argument requires Julia 1.4 or later.
Examples
julia> A = [3.0 -6.0; 4.0 -8.0; 0.0 1.0]3×2 Matrix{Float64}: 3.0 -6.0 4.0 -8.0 0.0 1.0julia> F = qr(A)LinearAlgebra.QRCompactWY{Float64, Matrix{Float64}, Matrix{Float64}}Q factor: 3×3 LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}, Matrix{Float64}}R factor:2×2 Matrix{Float64}: -5.0 10.0 0.0 -1.0julia> F.Q * F.R == Atrueqr returns multiple types because LAPACK uses several representations that minimize the memory storage requirements of products of Householder elementary reflectors, so that theQ andR matrices can be stored compactly rather as two separate dense matrices.
Several other Julia packages provide sparse matrix implementations that should be mentioned:
SuiteSparseGraphBLAS.jl is a wrapper over the fast, multithreaded SuiteSparse:GraphBLAS C library. On CPU this is typically the fastest option, often significantly outperforming MKLSparse.
CUDA.jl exposes theCUSPARSE library for GPU sparse matrix operations.
SparseMatricesCSR.jl provides a Julia native implementation of the Compressed Sparse Rows (CSR) format.
MKLSparse.jl accelerates SparseArrays sparse-dense matrix operations using Intel's MKL library.
SparseArrayKit.jl available for multidimensional sparse arrays.
LuxurySparse.jl provides static sparse array formats, as well as a coordinate format.
ExtendableSparse.jl enables fast insertion into sparse matrices using a lazy approach to new stored indices.
Finch.jl supports extensive multidimensional sparse array formats and operations through a mini tensor language and compiler, all in native Julia. Support for COO, CSF, CSR, CSC and more, as well as operations like broadcast, reduce, etc. and custom operations.
External packages providing sparse direct solvers:
External packages providing solvers for iterative solution of eigensystems and singular value decompositions:
External packages for working with graphs:
Settings
This document was generated withDocumenter.jl version 1.8.0 onWednesday 9 July 2025. Using Julia version 1.11.6.