Linear Algebra¶
Standard Functions¶
Linear algebra functions in Julia are largely implemented by calling functions from LAPACK. Sparse factorizations call functions from SuiteSparse.
-
*
(A, B)¶ Matrix multiplication
-
\
(A, B)¶ Matrix division using a polyalgorithm. For input matrices
A
andB
, the resultX
is such thatA*X == B
whenA
is square. The solver that is used depends upon the structure ofA
. A direct solver is used for upper or lower triangularA
. For HermitianA
(equivalent to symmetricA
for non-complexA
) theBunchKaufman
factorization is used. Otherwise an LU factorization is used. For rectangularA
the result is the minimum-norm least squares solution computed by a pivoted QR factorization ofA
and a rank estimate ofA
based on the R factor.When
A
is sparse, a similar polyalgorithm is used. For indefinite matrices, the LDLt factorization does not use pivoting during the numerical factorization and therefore the procedure can fail even for invertible matrices.
-
vecdot
(x, y)¶ For any iterable containers
x
andy
(including arrays of any dimension) of numbers (or any element type for whichdot
is defined), compute the Euclidean dot product (the sum ofdot(x[i],y[i])
) as if they were vectors.
-
factorize
(A)¶ Compute a convenient factorization (including LU, Cholesky, Bunch-Kaufman, LowerTriangular, UpperTriangular) of
A
, based upon the type of the input matrix. The return value can then be reused for efficient solving of multiple systems. For example:A=factorize(A); x=A\b; y=A\C
.
-
full
(F)¶ Reconstruct the matrix
A
from the factorizationF=factorize(A)
.
-
lu
(A) → L, U, p¶ Compute the LU factorization of
A
, such thatA[p,:] = L*U
.
-
lufact
(A[, pivot=Val{true}]) → F¶ Compute the LU factorization of
A
. The return type ofF
depends on the type ofA
. In most cases, ifA
is a subtypeS
of AbstractMatrix with an element typeT
supporting+
,-
,*
and/
the return type isLU{T,S{T}}
. If pivoting is chosen (default) the element type should also supportabs
and<
. WhenA
is sparse and have element of typeFloat32
,Float64
,Complex{Float32}
, orComplex{Float64}
the return type isUmfpackLU
. Some examples are shown in the table below.Type of input A
Type of output F
Relationship between F
andA
Matrix()
LU
F[:L]*F[:U] == A[F[:p], :]
Tridiagonal()
LU{T,Tridiagonal{T}}
F[:L]*F[:U] == A[F[:p], :]
SparseMatrixCSC()
UmfpackLU
F[:L]*F[:U] == (F[:Rs] .* A)[F[:p], F[:q]]
The individual components of the factorization
F
can be accessed by indexing:Component Description LU
LU{T,Tridiagonal{T}}
UmfpackLU
F[:L]
L
(lower triangular) part ofLU
✓ ✓ ✓ F[:U]
U
(upper triangular) part ofLU
✓ ✓ ✓ F[:p]
(right) permutation Vector
✓ ✓ ✓ F[:P]
(right) permutation Matrix
✓ ✓ F[:q]
left permutation Vector
✓ F[:Rs]
Vector
of scaling factors✓ F[:(:)]
(L,U,p,q,Rs)
components✓ Supported function LU
LU{T,Tridiagonal{T}}
UmfpackLU
/
✓ \
✓ ✓ ✓ cond
✓ ✓ det
✓ ✓ ✓ logdet
✓ ✓ logabsdet
✓ ✓ size
✓ ✓
-
lufact!
(A) → LU¶ lufact!
is the same aslufact()
, but saves space by overwriting the inputA
, instead of creating a copy. For sparseA
thenzval
field is not overwritten but the index fields,colptr
androwval
are decremented in place, converting from 1-based indices to 0-based indices.
-
chol
(A[, LU]) → F¶ Compute the Cholesky factorization of a symmetric positive definite matrix
A
and return the matrixF
. IfLU
isVal{:U}
(Upper),F
is of typeUpperTriangular
andA = F'*F
. IfLU
isVal{:L}
(Lower),F
is of typeLowerTriangular
andA = F*F'
.LU
defaults toVal{:U}
.
-
cholfact
(A, [LU=:U[,pivot=Val{false}]][;tol=-1.0]) → Cholesky¶ Compute the Cholesky factorization of a dense symmetric positive (semi)definite matrix
A
and return either aCholesky
ifpivot==Val{false}
orCholeskyPivoted
ifpivot==Val{true}
.LU
may be:L
for using the lower part or:U
for the upper part. The default is to use:U
. The triangular matrix can be obtained from the factorizationF
with:F[:L]
andF[:U]
. The following functions are available forCholesky
objects:size
,\
,inv
,det
. ForCholeskyPivoted
there is also defined arank
. Ifpivot==Val{false}
aPosDefException
exception is thrown in case the matrix is not positive definite. The argumenttol
determines the tolerance for determining the rank. For negative values, the tolerance is the machine precision.
-
cholfact
(A; shift=0, perm=Int[]) → CHOLMOD.Factor Compute the Cholesky factorization of a sparse positive definite matrix
A
. A fill-reducing permutation is used.F = cholfact(A)
is most frequently used to solve systems of equations withF\b
, but also the methodsdiag
,det
,logdet
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 extact “combined” factors likePtL = F[:PtL]
(the equivalent ofP'*L
) andLtP = F[:UP]
(the equivalent ofL'*P
).Setting optional
shift
keyword argument computes the factorization ofA+shift*I
instead ofA
. If theperm
argument is nonempty, it should be a permutation of1:size(A,1)
giving the ordering to use (instead of CHOLMOD’s default AMD ordering).The function calls the C library CHOLMOD and many other functions from the library are wrapped but not exported.
-
cholfact!
(A [,LU=:U [,pivot=Val{false}]][;tol=-1.0]) → Cholesky¶ cholfact!
is the same ascholfact()
, but saves space by overwriting the inputA
, instead of creating a copy.cholfact!
can also reuse the symbolic factorization from a different matrixF
with the same structure when used as:cholfact!(F::CholmodFactor, A)
.
-
ldltfact
(A) → LDLtFactorization¶ Compute a factorization of a positive definite matrix
A
such thatA=L*Diagonal(d)*L'
whereL
is a unit lower triangular matrix andd
is a vector with non-negative elements.
-
ldltfact
(A; shift=0, perm=Int[]) → CHOLMOD.Factor Compute the LDLt factorization of a sparse symmetric or Hermitian matrix
A
. A fill-reducing permutation is used.F = ldltfact(A)
is most frequently used to solve systems of equations withF\b
, but also the methodsdiag
,det
,logdet
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*D*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 extact “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
.Setting optional
shift
keyword argument computes the factorization ofA+shift*I
instead ofA
. If theperm
argument is nonempty, it should be a permutation of1:size(A,1)
giving the ordering to use (instead of CHOLMOD’s default AMD ordering).The function calls the C library CHOLMOD and many other functions from the library are wrapped but not exported.
-
qr
(A[, pivot=Val{false}][;thin=true]) → Q, R, [p]¶ Compute the (pivoted) QR factorization of
A
such that eitherA = Q*R
orA[:,p] = Q*R
. Also seeqrfact
. The default is to compute a thin factorization. Note thatR
is not extended with zeros when the fullQ
is requested.
-
qrfact
(A[, pivot=Val{false}]) → F¶ Computes the QR factorization of
A
. The return type ofF
depends on the element type ofA
and whether pivoting is specified (withpivot==Val{true}
).Return type eltype(A)
pivot
Relationship between F
andA
QR
not BlasFloat
either A==F[:Q]*F[:R]
QRCompactWY
BlasFloat
Val{false}
A==F[:Q]*F[:R]
QRPivoted
BlasFloat
Val{true}
A[:,F[:p]]==F[:Q]*F[:R]
BlasFloat
refers to any of:Float32
,Float64
,Complex64
orComplex128
.The individual components of the factorization
F
can be accessed by indexing:Component Description QR
QRCompactWY
QRPivoted
F[:Q]
Q
(orthogonal/unitary) part ofQR
✓ ( QRPackedQ
)✓ ( QRCompactWYQ
)✓ ( QRPackedQ
)F[:R]
R
(upper right triangular) part ofQR
✓ ✓ ✓ F[:p]
pivot Vector
✓ F[:P]
(pivot) permutation Matrix
✓ The following functions are available for the
QR
objects:size
,\
. WhenA
is rectangular,\
will return a least squares solution and if the solution is not unique, the one with smallest norm is returned.Multiplication with respect to either thin or full
Q
is allowed, i.e. bothF[:Q]*F[:R]
andF[:Q]*A
are supported. AQ
matrix can be converted into a regular matrix withfull()
which has a named argumentthin
.Note
qrfact
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.The data contained in
QR
orQRPivoted
can be used to construct theQRPackedQ
type, which is a compact representation of the rotation matrix:\[Q = \prod_{i=1}^{\min(m,n)} (I - \tau_i v_i v_i^T)\]where \(\tau_i\) is the scale factor and \(v_i\) is the projection vector associated with the \(i^{th}\) Householder elementary reflector.
The data contained in
QRCompactWY
can be used to construct theQRCompactWYQ
type, which is a compact representation of the rotation matrix\[Q = I + Y T Y^T\]where
Y
is \(m \times r\) lower trapezoidal andT
is \(r \times r\) upper triangular. The compact WY representation [Schreiber1989] is not to be confused with the older, WY representation [Bischof1987]. (The LAPACK documentation usesV
in lieu ofY
.)[Bischof1987] (1, 2) C Bischof and C Van Loan, “The WY representation for products of Householder matrices”, SIAM J Sci Stat Comput 8 (1987), s2-s13. doi:10.1137/0908009 [Schreiber1989] R Schreiber and C Van Loan, “A storage-efficient WY representation for products of Householder transformations”, SIAM J Sci Stat Comput 10 (1989), 53-57. doi:10.1137/0910005
-
qrfact
(A) → SPQR.Factorization Compute the QR factorization of a sparse matrix
A
. A fill-reducing permutation is used. The main application of this type is to solve least squares problems with\
. The function calls the C library SPQR and a few additional functions from the library are wrapped but not exported.
-
qrfact!
(A[, pivot=Val{false}])¶ qrfact!
is the same asqrfact()
whenA
is a subtype ofStridedMatrix
, but saves space by overwriting the inputA
, instead of creating a copy.
-
full
(QRCompactWYQ[, thin=true]) → Matrix Converts an orthogonal or unitary matrix stored as a
QRCompactWYQ
object, i.e. in the compact WY format [Bischof1987], to a dense matrix.Optionally takes a
thin
Boolean argument, which iftrue
omits the columns that span the rows ofR
in the QR factorization that are zero. The resulting matrix is theQ
in a thin QR factorization (sometimes called the reduced QR factorization). Iffalse
, returns aQ
that spans all rows ofR
in its corresponding QR factorization.
-
bkfact
(A) → BunchKaufman¶ Compute the Bunch-Kaufman [Bunch1977] factorization of a real symmetric or complex Hermitian matrix
A
and return aBunchKaufman
object. The following functions are available forBunchKaufman
objects:size
,\
,inv
,issym
,ishermitian
.[Bunch1977] J R Bunch and L Kaufman, Some stable methods for calculating inertia and solving symmetric linear systems, Mathematics of Computation 31:137 (1977), 163-179. url.
-
bkfact!
(A) → BunchKaufman¶ bkfact!
is the same asbkfact()
, but saves space by overwriting the inputA
, instead of creating a copy.
-
eig
(A,[irange,][vl,][vu,][permute=true,][scale=true]) → D, V¶ Computes eigenvalues and eigenvectors of
A
. Seeeigfact()
for details on thebalance
keyword argument.julia> eig([1.0 0.0 0.0; 0.0 3.0 0.0; 0.0 0.0 18.0]) ([1.0,3.0,18.0], 3x3 Array{Float64,2}: 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0)
eig
is a wrapper aroundeigfact()
, extracting all parts of the factorization to a tuple; where possible, usingeigfact()
is recommended.
-
eig
(A, B) → D, V Computes generalized eigenvalues and vectors of
A
with respect toB
.eig
is a wrapper aroundeigfact()
, extracting all parts of the factorization to a tuple; where possible, usingeigfact()
is recommended.
-
eigvals
(A,[irange,][vl,][vu])¶ Returns the eigenvalues of
A
. IfA
isSymmetric
,Hermitian
orSymTridiagonal
, it is possible to calculate only a subset of the eigenvalues by specifying either aUnitRange
irange
covering indices of the sorted eigenvalues, or a pairvl
andvu
for the lower and upper boundaries of the eigenvalues.For general non-symmetric matrices it is possible to specify how the matrix is balanced before the eigenvector calculation. The option
permute=true
permutes the matrix to become closer to upper triangular, andscale=true
scales the matrix by its diagonal elements to make rows and columns more equal in norm. The default istrue
for both options.
-
eigmax
(A)¶ Returns the largest eigenvalue of
A
.
-
eigmin
(A)¶ Returns the smallest eigenvalue of
A
.
-
eigvecs
(A, [eigvals,][permute=true,][scale=true]) → Matrix¶ Returns a matrix
M
whose columns are the eigenvectors ofA
. (Thek
th eigenvector can be obtained from the sliceM[:, k]
.) Thepermute
andscale
keywords are the same as foreigfact()
.For
SymTridiagonal
matrices, if the optional vector of eigenvalueseigvals
is specified, returns the specific corresponding eigenvectors.
-
eigfact
(A,[irange,][vl,][vu,][permute=true,][scale=true]) → Eigen¶ Computes the eigenvalue decomposition of
A
, returning anEigen
factorization objectF
which contains the eigenvalues inF[:values]
and the eigenvectors in the columns of the matrixF[:vectors]
. (Thek
th eigenvector can be obtained from the sliceF[:vectors][:, k]
.)The following functions are available for
Eigen
objects:inv
,det
.If
A
isSymmetric
,Hermitian
orSymTridiagonal
, it is possible to calculate only a subset of the eigenvalues by specifying either aUnitRange
irange
covering indices of the sorted eigenvalues or a pairvl
andvu
for the lower and upper boundaries of the eigenvalues.For general nonsymmetric matrices it is possible to specify how the matrix is balanced before the eigenvector calculation. The option
permute=true
permutes the matrix to become closer to upper triangular, andscale=true
scales the matrix by its diagonal elements to make rows and columns more equal in norm. The default istrue
for both options.
-
eigfact
(A, B) → GeneralizedEigen Computes the generalized eigenvalue decomposition of
A
andB
, returning aGeneralizedEigen
factorization objectF
which contains the generalized eigenvalues inF[:values]
and the generalized eigenvectors in the columns of the matrixF[:vectors]
. (Thek
th generalized eigenvector can be obtained from the sliceF[:vectors][:, k]
.)
-
eigfact!
(A[, B])¶ Same as
eigfact()
, but saves space by overwriting the inputA
(andB
), instead of creating a copy.
-
hessfact
(A)¶ Compute the Hessenberg decomposition of
A
and return aHessenberg
object. IfF
is the factorization object, the unitary matrix can be accessed withF[:Q]
and the Hessenberg matrix withF[:H]
. WhenQ
is extracted, the resulting type is theHessenbergQ
object, and may be converted to a regular matrix withfull()
.
-
hessfact!
(A)¶ hessfact!
is the same ashessfact()
, but saves space by overwriting the inputA
, instead of creating a copy.
-
schurfact
(A) → Schur¶ Computes the Schur factorization of the matrix
A
. The (quasi) triangular Schur factor can be obtained from theSchur
objectF
with eitherF[:Schur]
orF[:T]
and the unitary/orthogonal Schur vectors can be obtained withF[:vectors]
orF[:Z]
such thatA=F[:vectors]*F[:Schur]*F[:vectors]'
. The eigenvalues ofA
can be obtained withF[:values]
.
-
schurfact!
(A)¶ Computes the Schur factorization of
A
, overwritingA
in the process. Seeschurfact()
-
schur
(A) → Schur[:T], Schur[:Z], Schur[:values]¶ See
schurfact()
-
ordschur
(Q, T, select) → Schur¶ Reorders the Schur factorization of a real matrix
A=Q*T*Q'
according to the logical arrayselect
returning a Schur objectF
. The selected eigenvalues appear in the leading diagonal ofF[:Schur]
and the the corresponding leading columns ofF[:vectors]
form an orthonormal basis of the corresponding right invariant subspace. A complex conjugate pair of eigenvalues must be either both included or excluded viaselect
.
-
ordschur!
(Q, T, select) → Schur¶ Reorders the Schur factorization of a real matrix
A=Q*T*Q'
, overwritingQ
andT
in the process. Seeordschur()
-
ordschur
(S, select) → Schur Reorders the Schur factorization
S
of typeSchur
.
-
ordschur!
(S, select) → Schur Reorders the Schur factorization
S
of typeSchur
, overwritingS
in the process. Seeordschur()
-
schurfact
(A, B) → GeneralizedSchur Computes the Generalized Schur (or QZ) factorization of the matrices
A
andB
. The (quasi) triangular Schur factors can be obtained from theSchur
objectF
withF[:S]
andF[:T]
, the left unitary/orthogonal Schur vectors can be obtained withF[:left]
orF[:Q]
and the right unitary/orthogonal Schur vectors can be obtained withF[:right]
orF[:Z]
such thatA=F[:left]*F[:S]*F[:right]'
andB=F[:left]*F[:T]*F[:right]'
. The generalized eigenvalues ofA
andB
can be obtained withF[:alpha]./F[:beta]
.
-
schur
(A, B) → GeneralizedSchur[:S], GeneralizedSchur[:T], GeneralizedSchur[:Q], GeneralizedSchur[:Z] See
schurfact()
-
ordschur
(S, T, Q, Z, select) → GeneralizedSchur Reorders the Generalized Schur factorization of a matrix
(A, B) = (Q*S*Z^{H}, Q*T*Z^{H})
according to the logical arrayselect
and returns a GeneralizedSchur objectGS
. The selected eigenvalues appear in the leading diagonal of both(GS[:S], GS[:T])
and the left and right unitary/orthogonal Schur vectors are also reordered such that(A, B) = GS[:Q]*(GS[:S], GS[:T])*GS[:Z]^{H}
still holds and the generalized eigenvalues ofA
andB
can still be obtained withGS[:alpha]./GS[:beta]
.
-
ordschur!
(S, T, Q, Z, select) → GeneralizedSchur Reorders the Generalized Schur factorization of a matrix by overwriting the matrices
(S, T, Q, Z)
in the process. Seeordschur()
.
-
ordschur
(GS, select) → GeneralizedSchur Reorders the Generalized Schur factorization of a Generalized Schur object. See
ordschur()
.
-
ordschur!
(GS, select) → GeneralizedSchur Reorders the Generalized Schur factorization of a Generalized Schur object by overwriting the object with the new factorization. See
ordschur()
.
-
svdfact
(A[, thin=true]) → SVD¶ Compute the Singular Value Decomposition (SVD) of
A
and return anSVD
object.U
,S
,V
andVt
can be obtained from the factorizationF
withF[:U]
,F[:S]
,F[:V]
andF[:Vt]
, such thatA = U*diagm(S)*Vt
. Ifthin
istrue
, an economy mode decomposition is returned. The algorithm producesVt
and henceVt
is more efficient to extract thanV
. The default is to produce a thin decomposition.
-
svdfact!
(A[, thin=true]) → SVD¶ svdfact!
is the same assvdfact()
, but saves space by overwriting the inputA
, instead of creating a copy. Ifthin
istrue
, an economy mode decomposition is returned. The default is to produce a thin decomposition.
-
svd
(A[, thin=true]) → U, S, V¶ Wrapper around
svdfact
extracting all parts the factorization to a tuple. Direct use ofsvdfact
is therefore generally more efficient. Computes the SVD ofA
, returningU
, vectorS
, andV
such thatA == U*diagm(S)*V'
. Ifthin
istrue
, an economy mode decomposition is returned. The default is to produce a thin decomposition.
-
svdvals
(A)¶ Returns the singular values of
A
.
-
svdvals!
(A)¶ Returns the singular values of
A
, while saving space by overwriting the input.
-
svdfact
(A, B) → GeneralizedSVD Compute the generalized SVD of
A
andB
, returning aGeneralizedSVD
Factorization objectF
, such thatA = F[:U]*F[:D1]*F[:R0]*F[:Q]'
andB = F[:V]*F[:D2]*F[:R0]*F[:Q]'
.
-
svd
(A, B) → U, V, Q, D1, D2, R0 Wrapper around
svdfact
extracting all parts the factorization to a tuple. Direct use ofsvdfact
is therefore generally more efficient. The function returns the generalized SVD ofA
andB
, returningU
,V
,Q
,D1
,D2
, andR0
such thatA = U*D1*R0*Q'
andB = V*D2*R0*Q'
.
-
svdvals
(A, B) Return only the singular values from the generalized singular value decomposition of
A
andB
.
-
triu
(M)¶ Upper triangle of a matrix.
-
triu
(M, k) Returns the upper triangle of
M
starting from thek
th superdiagonal.
-
triu!
(M)¶ Upper triangle of a matrix, overwriting
M
in the process.
-
triu!
(M, k) Returns the upper triangle of
M
starting from thek
th superdiagonal, overwritingM
in the process.
-
tril
(M)¶ Lower triangle of a matrix.
-
tril
(M, k) Returns the lower triangle of
M
starting from thek
th superdiagonal.
-
tril!
(M)¶ Lower triangle of a matrix, overwriting
M
in the process.
-
tril!
(M, k) Returns the lower triangle of
M
starting from thek
th superdiagonal, overwritingM
in the process.
-
diagind
(M[, k])¶ A
Range
giving the indices of thek
th diagonal of the matrixM
.
-
diag
(M[, k])¶ The
k
th diagonal of a matrix, as a vector. Usediagm
to construct a diagonal matrix.
-
diagm
(v[, k])¶ Construct a diagonal matrix and place
v
on thek
th diagonal.
-
scale
(A, b)¶ -
scale
(b, A) Scale an array
A
by a scalarb
, returning a new array.If
A
is a matrix andb
is a vector, thenscale(A,b)
scales each columni
ofA
byb[i]
(similar toA*diagm(b)
), whilescale(b,A)
scales each rowi
ofA
byb[i]
(similar todiagm(b)*A
), returning a new array.Note: for large
A
,scale
can be much faster thanA .* b
orb .* A
, due to the use of BLAS.
-
scale!
(A, b)¶ -
scale!
(b, A) Scale an array
A
by a scalarb
, similar toscale()
but overwritingA
in-place.If
A
is a matrix andb
is a vector, thenscale!(A,b)
scales each columni
ofA
byb[i]
(similar toA*diagm(b)
), whilescale!(b,A)
scales each rowi
ofA
byb[i]
(similar todiagm(b)*A
), again operating in-place onA
.
-
Tridiagonal
(dl, d, du)¶ Construct a tridiagonal matrix from the lower diagonal, diagonal, and upper diagonal, respectively. The result is of type
Tridiagonal
and provides efficient specialized linear solvers, but may be converted into a regular matrix withfull()
.
-
Bidiagonal
(dv, ev, isupper)¶ Constructs an upper (
isupper=true
) or lower (isupper=false
) bidiagonal matrix using the given diagonal (dv
) and off-diagonal (ev
) vectors. The result is of typeBidiagonal
and provides efficient specialized linear solvers, but may be converted into a regular matrix withfull()
.
-
SymTridiagonal
(d, du)¶ Construct a real symmetric tridiagonal matrix from the diagonal and upper diagonal, respectively. The result is of type
SymTridiagonal
and provides efficient specialized eigensolvers, but may be converted into a regular matrix withfull()
.
-
rank
(M)¶ Compute the rank of a matrix.
-
norm
(A[, p])¶ Compute the
p
-norm of a vector or the operator norm of a matrixA
, defaulting to thep=2
-norm.For vectors,
p
can assume any numeric value (even though not all values produce a mathematically valid vector norm). In particular,norm(A, Inf)
returns the largest value inabs(A)
, whereasnorm(A, -Inf)
returns the smallest.For matrices, the matrix norm induced by the vector
p
-norm is used, where valid values ofp
are1
,2
, orInf
. (Note that for sparse matrices,p=2
is currently not implemented.) Usevecnorm()
to compute the Frobenius norm.
-
vecnorm
(A[, p])¶ For any iterable container
A
(including arrays of any dimension) of numbers (or any element type for whichnorm
is defined), compute thep
-norm (defaulting top=2
) as ifA
were a vector of the corresponding length.For example, if
A
is a matrix andp=2
, then this is equivalent to the Frobenius norm.
-
cond
(M[, p])¶ Condition number of the matrix
M
, computed using the operatorp
-norm. Valid values forp
are1
,2
(default), orInf
.
-
condskeel
(M[, x, p])¶ - \[\begin{split}\kappa_S(M, p) & = \left\Vert \left\vert M \right\vert \left\vert M^{-1} \right\vert \right\Vert_p \\ \kappa_S(M, x, p) & = \left\Vert \left\vert M \right\vert \left\vert M^{-1} \right\vert \left\vert x \right\vert \right\Vert_p\end{split}\]
Skeel condition number \(\kappa_S\) of the matrix
M
, optionally with respect to the vectorx
, as computed using the operatorp
-norm.p
isInf
by default, if not provided. Valid values forp
are1
,2
, orInf
.This quantity is also known in the literature as the Bauer condition number, relative condition number, or componentwise relative condition number.
-
trace
(M)¶ Matrix trace
-
det
(M)¶ Matrix determinant
-
logdet
(M)¶ Log of matrix determinant. Equivalent to
log(det(M))
, but may provide increased accuracy and/or speed.
-
logabsdet
(M)¶ Log of absolute value of determinant of real matrix. Equivalent to
(log(abs(det(M))), sign(det(M)))
, but may provide increased accuracy and/or speed.
-
inv
(M)¶ Matrix inverse
-
pinv
(M[, tol])¶ Computes the Moore-Penrose pseudoinverse.
For matrices
M
with floating point elements, it is convenient to compute the pseudoinverse by inverting only singular values above a given threshold,tol
.The optimal choice of
tol
varies both with the value ofM
and the intended application of the pseudoinverse. The default value oftol
iseps(real(float(one(eltype(M)))))*maximum(size(A))
, which is essentially machine epsilon for the real part of a matrix element multiplied by the larger matrix dimension. For inverting dense ill-conditioned matrices in a least-squares sense,tol = sqrt(eps(real(float(one(eltype(M))))))
is recommended.For more information, see [8859], [B96], [S84], [KY88].
[8859] Issue 8859, “Fix least squares”, https://github.com/JuliaLang/julia/pull/8859 [B96] Åke Björck, “Numerical Methods for Least Squares Problems”, SIAM Press, Philadelphia, 1996, “Other Titles in Applied Mathematics”, Vol. 51. doi:10.1137/1.9781611971484 [S84] G. W. Stewart, “Rank Degeneracy”, SIAM Journal on Scientific and Statistical Computing, 5(2), 1984, 403-413. doi:10.1137/0905030 [KY88] Konstantinos Konstantinides and Kung Yao, “Statistical analysis of effective singular values in matrix rank determination”, IEEE Transactions on Acoustics, Speech and Signal Processing, 36(5), 1988, 757-763. doi:10.1109/29.1585
-
nullspace
(M)¶ Basis for nullspace of
M
.
-
repmat
(A, n, m)¶ Construct a matrix by repeating the given matrix
n
times in dimension 1 andm
times in dimension 2.
-
repeat
(A, inner = Int[], outer = Int[])¶ Construct an array by repeating the entries of
A
. The i-th element ofinner
specifies the number of times that the individual entries of the i-th dimension ofA
should be repeated. The i-th element ofouter
specifies the number of times that a slice along the i-th dimension ofA
should be repeated.
-
kron
(A, B)¶ Kronecker tensor product of two vectors or two matrices.
-
blkdiag
(A...)¶ Concatenate matrices block-diagonally. Currently only implemented for sparse matrices.
-
linreg
(x, y) → [a; b]¶ Linear Regression. Returns
a
andb
such thata+b*x
is the closest line to the given points(x,y)
. In other words, this function determines parameters[a, b]
that minimize the squared error betweeny
anda+b*x
.Example:
using PyPlot; x = float([1:12]) y = [5.5; 6.3; 7.6; 8.8; 10.9; 11.79; 13.48; 15.02; 17.77; 20.81; 22.0; 22.99] a, b = linreg(x,y) # Linear regression plot(x, y, "o") # Plot (x,y) points plot(x, [a+b*i for i in x]) # Plot the line determined by the linear regression
-
linreg
(x, y, w) Weighted least-squares linear regression.
-
expm
(A)¶ Compute the matrix exponential of
A
, defined by\[e^A = \sum_{n=0}^{\infty} \frac{A^n}{n!}.\]For symmetric or Hermitian
A
, an eigendecomposition (eigfact()
) is used, otherwise the scaling and squaring algorithm (see [H05]) is chosen.[H05] Nicholas J. Higham, “The squaring and scaling method for the matrix exponential revisited”, SIAM Journal on Matrix Analysis and Applications, 26(4), 2005, 1179-1193. doi:10.1137/090768539
-
logm
(A)¶ If
A
has no negative real eigenvalue, compute the principal matrix logarithm ofA
, i.e. the unique matrix \(X\) such that \(e^X = A\) and \(-\pi < Im(\lambda) < \pi\) for all the eigenvalues \(\lambda\) of \(X\). IfA
has nonpositive eigenvalues, a warning is printed and whenever possible a nonprincipal matrix function is returned.If
A
is symmetric or Hermitian, its eigendecomposition (eigfact()
) is used, ifA
is triangular an improved version of the inverse scaling and squaring method is employed (see [AH12] and [AHR13]). For general matrices, the complex Schur form (schur()
) is computed and the triangular algorithm is used on the triangular factor.[AH12] Awad H. Al-Mohy and Nicholas J. Higham, “Improved inverse scaling and squaring algorithms for the matrix logarithm”, SIAM Journal on Scientific Computing, 34(4), 2012, C153-C169. doi:10.1137/110852553 [AHR13] Awad H. Al-Mohy, Nicholas J. Higham and Samuel D. Relton, “Computing the Fréchet derivative of the matrix logarithm and estimating the condition number”, SIAM Journal on Scientific Computing, 35(4), 2013, C394-C410. doi:10.1137/120885991
-
sqrtm
(A)¶ If
A
has no negative real eigenvalues, compute the principal matrix square root ofA
, that is the unique matrix \(X\) with eigenvalues having positive real part such that \(X^2 = A\). Otherwise, a nonprincipal square root is returned.If
A
is symmetric or Hermitian, its eigendecomposition (eigfact()
) is used to compute the square root. Otherwise, the square root is determined by means of the Björck-Hammarling method, which computes the complex Schur form (schur()
) and then the complex square root of the triangular factor.[BH83] Åke Björck and Sven Hammarling, “A Schur method for the square root of a matrix”, Linear Algebra and its Applications, 52-53, 1983, 127-140. doi:10.1016/0024-3795(83)80010-X
-
lyap
(A, C)¶ Computes the solution
X
to the continuous Lyapunov equationAX + XA' + C = 0
, where no eigenvalue ofA
has a zero real part and no two eigenvalues are negative complex conjugates of each other.
-
sylvester
(A, B, C)¶ Computes the solution
X
to the Sylvester equationAX + XB + C = 0
, whereA
,B
andC
have compatible dimensions andA
and-B
have no eigenvalues with equal real part.
-
issym
(A) → Bool¶ Test whether a matrix is symmetric.
-
isposdef
(A) → Bool¶ Test whether a matrix is positive definite.
-
isposdef!
(A) → Bool¶ Test whether a matrix is positive definite, overwriting
A
in the processes.
-
istril
(A) → Bool¶ Test whether a matrix is lower triangular.
-
istriu
(A) → Bool¶ Test whether a matrix is upper triangular.
-
isdiag
(A) → Bool¶ Test whether a matrix is diagonal.
-
ishermitian
(A) → Bool¶ Test whether a matrix is Hermitian.
-
transpose
(A)¶ The transposition operator (
.'
).
-
transpose!
(dest, src)¶ Transpose array
src
and store the result in the preallocated arraydest
, which should have a size corresponding to(size(src,2),size(src,1))
. No in-place transposition is supported and unexpected results will happen ifsrc
anddest
have overlapping memory regions.
-
ctranspose
(A)¶ The conjugate transposition operator (
'
).
-
ctranspose!
(dest, src)¶ Conjugate transpose array
src
and store the result in the preallocated arraydest
, which should have a size corresponding to(size(src,2),size(src,1))
. No in-place transposition is supported and unexpected results will happen ifsrc
anddest
have overlapping memory regions.
-
eigs
(A, [B, ]; nev=6, which="LM", tol=0.0, maxiter=300, sigma=nothing, ritzvec=true, v0=zeros((0, ))) -> (d, [v, ]nconv, niter, nmult, resid)¶ Computes eigenvalues
d
ofA
using Lanczos or Arnoldi iterations for real symmetric or general nonsymmetric matrices respectively. IfB
is provided, the generalized eigenproblem is solved.- The following keyword arguments are supported:
nev
: Number of eigenvaluesncv
: Number of Krylov vectors used in the computation; should satisfynev+1 <= ncv <= n
for real symmetric problems andnev+2 <= ncv <= n
for other problems, wheren
is the size of the input matrixA
. The default isncv = max(20,2*nev+1)
.Note that these restrictions limit the input matrix
A
to be of dimension at least 2.which
: type of eigenvalues to compute. See the note below.which
type of eigenvalues :LM
eigenvalues of largest magnitude (default) :SM
eigenvalues of smallest magnitude :LR
eigenvalues of largest real part :SR
eigenvalues of smallest real part :LI
eigenvalues of largest imaginary part (nonsymmetric or complex A
only):SI
eigenvalues of smallest imaginary part (nonsymmetric or complex A
only):BE
compute half of the eigenvalues from each end of the spectrum, biased in favor of the high end. (real symmetric A
only)tol
: tolerance (\(tol \le 0.0\) defaults toDLAMCH('EPS')
)maxiter
: Maximum number of iterations (default = 300)sigma
: Specifies the level shift used in inverse iteration. Ifnothing
(default), defaults to ordinary (forward) iterations. Otherwise, find eigenvalues close tosigma
using shift and invert iterations.ritzvec
: Returns the Ritz vectorsv
(eigenvectors) iftrue
v0
: starting vector from which to start the iterations
eigs
returns thenev
requested eigenvalues ind
, the corresponding Ritz vectorsv
(only ifritzvec=true
), the number of converged eigenvaluesnconv
, the number of iterationsniter
and the number of matrix vector multiplicationsnmult
, as well as the final residual vectorresid
.Note
The
sigma
andwhich
keywords interact: the description of eigenvalues searched for bywhich
do _not_ necessarily refer to the eigenvalues ofA
, but rather the linear operator constructed by the specification of the iteration mode implied bysigma
.sigma
iteration mode which
refers to eigenvalues ofnothing
ordinary (forward) \(A\) real or complex inverse with level shift sigma
\((A - \sigma I )^{-1}\)
-
svds
(A; nsv=6, ritzvec=true, tol=0.0, maxiter=1000) -> (left_sv, s, right_sv, nconv, niter, nmult, resid)¶ svds
computes largest singular valuess
ofA
using Lanczos or Arnoldi iterations. Useseigs()
underneath.- Inputs are:
A
: Linear operator. It can either subtype ofAbstractArray
(e.g., sparse matrix) or duck typed. For duck typingA
has to supportsize(A)
,eltype(A)
,A * vector
andA' * vector
.nsv
: Number of singular values.ritzvec
: Whether to return the left and right singular vectorsleft_sv
andright_sv
, default istrue
. Iffalse
the singular vectors are omitted from the output.tol
: tolerance, seeeigs()
.maxiter
: Maximum number of iterations, seeeigs()
.
Example:
X = sprand(10, 5, 0.2) svds(X, nsv = 2)
-
peakflops
(n; parallel=false)¶ peakflops
computes the peak flop rate of the computer by using double precisionBase.LinAlg.BLAS.gemm!()
. By default, if no arguments are specified, it multiplies a matrix of sizen x n
, wheren = 2000
. If the underlying BLAS is using multiple threads, higher flop rates are realized. The number of BLAS threads can be set withblas_set_num_threads(n)
.If the keyword argument
parallel
is set totrue
,peakflops
is run in parallel on all the worker processors. The flop rate of the entire parallel computer is returned. When running in parallel, only 1 BLAS thread is used. The argumentn
still refers to the size of the problem that is solved on each processor.
BLAS Functions¶
Base.LinAlg.BLAS
provides wrappers for some of the BLAS functions for
linear algebra. Those BLAS functions that overwrite one of the input
arrays have names ending in '!'
.
Usually a function has 4 methods defined, one each for Float64
,
Float32
, Complex128
and Complex64
arrays.
-
dot
(n, X, incx, Y, incy)¶ Dot product of two vectors consisting of
n
elements of arrayX
with strideincx
andn
elements of arrayY
with strideincy
.
-
dotu
(n, X, incx, Y, incy)¶ Dot function for two complex vectors.
-
dotc
(n, X, incx, U, incy)¶ Dot function for two complex vectors conjugating the first vector.
-
blascopy!
(n, X, incx, Y, incy)¶ Copy
n
elements of arrayX
with strideincx
to arrayY
with strideincy
. ReturnsY
.
-
nrm2
(n, X, incx)¶ 2-norm of a vector consisting of
n
elements of arrayX
with strideincx
.
-
asum
(n, X, incx)¶ sum of the absolute values of the first
n
elements of arrayX
with strideincx
.
-
axpy!
(a, X, Y)¶ Overwrite
Y
witha*X + Y
. ReturnsY
.
-
scal!
(n, a, X, incx)¶ Overwrite
X
witha*X
. ReturnsX
.
-
scal
(n, a, X, incx)¶ Returns
a*X
.
-
ger!
(alpha, x, y, A)¶ Rank-1 update of the matrix
A
with vectorsx
andy
asalpha*x*y' + A
.
-
syr!
(uplo, alpha, x, A)¶ Rank-1 update of the symmetric matrix
A
with vectorx
asalpha*x*x.' + A
. Whenuplo
is ‘U’ the upper triangle ofA
is updated (‘L’ for lower triangle). ReturnsA
.
-
syrk!
(uplo, trans, alpha, A, beta, C)¶ Rank-k update of the symmetric matrix
C
asalpha*A*A.' + beta*C
oralpha*A.'*A + beta*C
according to whethertrans
is ‘N’ or ‘T’. Whenuplo
is ‘U’ the upper triangle ofC
is updated (‘L’ for lower triangle). ReturnsC
.
-
syrk
(uplo, trans, alpha, A)¶ Returns either the upper triangle or the lower triangle, according to
uplo
(‘U’ or ‘L’), ofalpha*A*A.'
oralpha*A.'*A
, according totrans
(‘N’ or ‘T’).
-
her!
(uplo, alpha, x, A)¶ Methods for complex arrays only. Rank-1 update of the Hermitian matrix
A
with vectorx
asalpha*x*x' + A
. Whenuplo
is ‘U’ the upper triangle ofA
is updated (‘L’ for lower triangle). ReturnsA
.
-
herk!
(uplo, trans, alpha, A, beta, C)¶ Methods for complex arrays only. Rank-k update of the Hermitian matrix
C
asalpha*A*A' + beta*C
oralpha*A'*A + beta*C
according to whethertrans
is ‘N’ or ‘T’. Whenuplo
is ‘U’ the upper triangle ofC
is updated (‘L’ for lower triangle). ReturnsC
.
-
herk
(uplo, trans, alpha, A)¶ Methods for complex arrays only. Returns either the upper triangle or the lower triangle, according to
uplo
(‘U’ or ‘L’), ofalpha*A*A'
oralpha*A'*A
, according totrans
(‘N’ or ‘T’).
-
gbmv!
(trans, m, kl, ku, alpha, A, x, beta, y)¶ Update vector
y
asalpha*A*x + beta*y
oralpha*A'*x + beta*y
according totrans
(‘N’ or ‘T’). The matrixA
is a general band matrix of dimensionm
bysize(A,2)
withkl
sub-diagonals andku
super-diagonals. Returns the updatedy
.
-
gbmv
(trans, m, kl, ku, alpha, A, x, beta, y)¶ Returns
alpha*A*x
oralpha*A'*x
according totrans
(‘N’ or ‘T’). The matrixA
is a general band matrix of dimensionm
bysize(A,2)
withkl
sub-diagonals andku
super-diagonals.
-
sbmv!
(uplo, k, alpha, A, x, beta, y)¶ Update vector
y
asalpha*A*x + beta*y
whereA
is a a symmetric band matrix of ordersize(A,2)
withk
super-diagonals stored in the argumentA
. The storage layout forA
is described the reference BLAS module, level-2 BLAS at <http://www.netlib.org/lapack/explore-html/>.Returns the updated
y
.
-
sbmv
(uplo, k, alpha, A, x)¶ Returns
alpha*A*x
whereA
is a symmetric band matrix of ordersize(A,2)
withk
super-diagonals stored in the argumentA
.
-
sbmv
(uplo, k, A, x) Returns
A*x
whereA
is a symmetric band matrix of ordersize(A,2)
withk
super-diagonals stored in the argumentA
.
-
gemm!
(tA, tB, alpha, A, B, beta, C)¶ Update
C
asalpha*A*B + beta*C
or the other three variants according totA
(transposeA
) andtB
. Returns the updatedC
.
-
gemm
(tA, tB, alpha, A, B)¶ Returns
alpha*A*B
or the other three variants according totA
(transposeA
) andtB
.
-
gemm
(tA, tB, A, B) Returns
A*B
or the other three variants according totA
(transposeA
) andtB
.
-
gemv!
(tA, alpha, A, x, beta, y)¶ Update the vector
y
asalpha*A*x + beta*y
oralpha*A'x + beta*y
according totA
(transposeA
). Returns the updatedy
.
-
gemv
(tA, alpha, A, x)¶ Returns
alpha*A*x
oralpha*A'x
according totA
(transposeA
).
-
gemv
(tA, A, x) Returns
A*x
orA'x
according totA
(transposeA
).
-
symm!
(side, ul, alpha, A, B, beta, C)¶ Update
C
asalpha*A*B + beta*C
oralpha*B*A + beta*C
according toside
.A
is assumed to be symmetric. Only theul
triangle ofA
is used. Returns the updatedC
.
-
symm
(side, ul, alpha, A, B)¶ Returns
alpha*A*B
oralpha*B*A
according toside
.A
is assumed to be symmetric. Only theul
triangle ofA
is used.
-
symm
(side, ul, A, B) Returns
A*B
orB*A
according toside
.A
is assumed to be symmetric. Only theul
triangle ofA
is used.
-
symm
(tA, tB, alpha, A, B) Returns
alpha*A*B
or the other three variants according totA
(transposeA
) andtB
.
-
symv!
(ul, alpha, A, x, beta, y)¶ Update the vector
y
asalpha*A*x + beta*y
.A
is assumed to be symmetric. Only theul
triangle ofA
is used. Returns the updatedy
.
-
symv
(ul, alpha, A, x)¶ Returns
alpha*A*x
.A
is assumed to be symmetric. Only theul
triangle ofA
is used.
-
symv
(ul, A, x) Returns
A*x
.A
is assumed to be symmetric. Only theul
triangle ofA
is used.
-
trmm!
(side, ul, tA, dA, alpha, A, B)¶ Update
B
asalpha*A*B
or one of the other three variants determined byside
(A
on left or right) andtA
(transposeA
). Only theul
triangle ofA
is used.dA
indicates ifA
is unit-triangular (the diagonal is assumed to be all ones). Returns the updatedB
.
-
trmm
(side, ul, tA, dA, alpha, A, B)¶ Returns
alpha*A*B
or one of the other three variants determined byside
(A
on left or right) andtA
(transposeA
). Only theul
triangle ofA
is used.dA
indicates ifA
is unit-triangular (the diagonal is assumed to be all ones).
-
trsm!
(side, ul, tA, dA, alpha, A, B)¶ Overwrite
B
with the solution toA*X = alpha*B
or one of the other three variants determined byside
(A
on left or right ofX
) andtA
(transposeA
). Only theul
triangle ofA
is used.dA
indicates ifA
is unit-triangular (the diagonal is assumed to be all ones). Returns the updatedB
.
-
trsm
(side, ul, tA, dA, alpha, A, B)¶ Returns the solution to
A*X = alpha*B
or one of the other three variants determined byside
(A
on left or right ofX
) andtA
(transposeA
). Only theul
triangle ofA
is used.dA
indicates ifA
is unit-triangular (the diagonal is assumed to be all ones).
-
trmv!
(side, ul, tA, dA, alpha, A, b)¶ Update
b
asalpha*A*b
or one of the other three variants determined byside
(A
on left or right) andtA
(transposeA
). Only theul
triangle ofA
is used.dA
indicates ifA
is unit-triangular (the diagonal is assumed to be all ones). Returns the updatedb
.
-
trmv
(side, ul, tA, dA, alpha, A, b)¶ Returns
alpha*A*b
or one of the other three variants determined byside
(A
on left or right) andtA
(transposeA
). Only theul
triangle ofA
is used.dA
indicates ifA
is unit-triangular (the diagonal is assumed to be all ones).
-
trsv!
(ul, tA, dA, A, b)¶ Overwrite
b
with the solution toA*x = b
or one of the other two variants determined bytA
(transposeA
) andul
(triangle ofA
used).dA
indicates ifA
is unit-triangular (the diagonal is assumed to be all ones). Returns the updatedb
.
-
trsv
(ul, tA, dA, A, b)¶ Returns the solution to
A*x = b
or one of the other two variants determined bytA
(transposeA
) andul
(triangle ofA
is used.)dA
indicates ifA
is unit-triangular (the diagonal is assumed to be all ones).
-
blas_set_num_threads
(n)¶ Set the number of threads the BLAS library should use.
-
I
¶ An object of type
UniformScaling
, representing an identity matrix of any size.