Linear Algebra

In addition to (and as part of) its support for multi-dimensional arrays, Julia provides native implementations of many common and useful linear algebra operations which can be loaded with using LinearAlgebra. Basic operations, such as tr, det, and inv are all supported:

julia> A = [1 2 3; 4 1 6; 7 8 1]
3×3 Array{Int64,2}:
 1  2  3
 4  1  6
 7  8  1
julia> tr(A)
3
julia> det(A)
104.0
julia> inv(A)
3×3 Array{Float64,2}:
 -0.451923   0.211538    0.0865385
  0.365385  -0.192308    0.0576923
  0.240385   0.0576923  -0.0673077

As well as other useful operations, such as finding eigenvalues or eigenvectors:

julia> A = [-4. -17.; 2. 2.]
2×2 Array{Float64,2}:
 -4.0  -17.0
  2.0    2.0
julia> eigvals(A)
2-element Array{Complex{Float64},1}:
 -1.0 + 5.0im
 -1.0 - 5.0im
julia> eigvecs(A)
2×2 Array{Complex{Float64},2}:
  0.945905+0.0im        0.945905-0.0im
 -0.166924-0.278207im  -0.166924+0.278207im

In addition, Julia provides many factorizations which can be used to speed up problems such as linear solve or matrix exponentiation by pre-factorizing a matrix into a form more amenable (for performance or memory reasons) to the problem. See the documentation on factorize for more information. As an example:

julia> A = [1.5 2 -4; 3 -1 -6; -10 2.3 4]
3×3 Array{Float64,2}:
   1.5   2.0  -4.0
   3.0  -1.0  -6.0
 -10.0   2.3   4.0
julia> factorize(A)
LU{Float64,Array{Float64,2}}
L factor:
3×3 Array{Float64,2}:
  1.0    0.0       0.0
 -0.15   1.0       0.0
 -0.3   -0.132196  1.0
U factor:
3×3 Array{Float64,2}:
 -10.0  2.3     4.0
   0.0  2.345  -3.4
   0.0  0.0    -5.24947

Since A is not Hermitian, symmetric, triangular, tridiagonal, or bidiagonal, an LU factorization may be the best we can do. Compare with:

julia> B = [1.5 2 -4; 2 -1 -3; -4 -3 5]
3×3 Array{Float64,2}:
  1.5   2.0  -4.0
  2.0  -1.0  -3.0
 -4.0  -3.0   5.0
julia> factorize(B)
BunchKaufman{Float64,Array{Float64,2}}
D factor:
3×3 Tridiagonal{Float64,Array{Float64,1}}:
 -1.64286   0.0   ⋅
  0.0      -2.8  0.0
   ⋅        0.0  5.0
U factor:
3×3 UnitUpperTriangular{Float64,Array{Float64,2}}:
 1.0  0.142857  -0.8
  ⋅   1.0       -0.6
  ⋅    ⋅         1.0
permutation:
3-element Array{Int64,1}:
 1
 2
 3

Here, Julia was able to detect that B is in fact symmetric, and used a more appropriate factorization. Often it’s possible to write more efficient code for a matrix that is known to have certain properties e.g. it is symmetric, or tridiagonal. Julia provides some special types so that you can “tag” matrices as having these properties. For instance:

julia> B = [1.5 2 -4; 2 -1 -3; -4 -3 5]
3×3 Array{Float64,2}:
  1.5   2.0  -4.0
  2.0  -1.0  -3.0
 -4.0  -3.0   5.0
julia> sB = Symmetric(B)
3×3 Symmetric{Float64,Array{Float64,2}}:
  1.5   2.0  -4.0
  2.0  -1.0  -3.0
 -4.0  -3.0   5.0

sB has been tagged as a matrix that’s (real) symmetric, so for later operations we might perform on it, such as eigenfactorization or computing matrix-vector products, efficiencies can be found by only referencing half of it. For example:

julia> B = [1.5 2 -4; 2 -1 -3; -4 -3 5]
3×3 Array{Float64,2}:
  1.5   2.0  -4.0
  2.0  -1.0  -3.0
 -4.0  -3.0   5.0
julia> sB = Symmetric(B)
3×3 Symmetric{Float64,Array{Float64,2}}:
  1.5   2.0  -4.0
  2.0  -1.0  -3.0
 -4.0  -3.0   5.0
julia> x = [1; 2; 3]
3-element Array{Int64,1}:
 1
 2
 3
julia> sB\x
3-element Array{Float64,1}:
 -1.7391304347826084
 -1.1086956521739126
 -1.4565217391304346

The \ operation here performs the linear solution. The left-division operator is pretty powerful and it’s easy to write compact, readable code that is flexible enough to solve all sorts of systems of linear equations.

Special matrices

Matrices with special symmetries and structures arise often in linear algebra and are frequently associated with various matrix factorizations. Julia features a rich collection of special matrix types, which allow for fast computation with specialized routines that are specially developed for particular matrix types.

The following tables summarize the types of special matrices that have been implemented in Julia, as well as whether hooks to various optimized methods for them in LAPACK are available.

Elementary operations

Legend:

Matrix factorizations

Legend:

The uniform scaling operator

A UniformScaling operator represents a scalar times the identity operator, λ*I. The identity operator I is defined as a constant and is an instance of UniformScaling. The size of these operators are generic and match the other matrix in the binary operations +, -, * and \. For A+I and A-I this means that A must be square. Multiplication with the identity operator I is a noop (except for checking that the scaling factor is one) and therefore almost without overhead.

To see the UniformScaling operator in action:

julia> U = UniformScaling(2);
julia> a = [1 2; 3 4]
2×2 Array{Int64,2}:
 1  2
 3  4
julia> a + U
2×2 Array{Int64,2}:
 3  2
 3  6
julia> a * U
2×2 Array{Int64,2}:
 2  4
 6  8
julia> [a U]
2×4 Array{Int64,2}:
 1  2  2  0
 3  4  0  2
julia> b = [1 2 3; 4 5 6]
2×3 Array{Int64,2}:
 1  2  3
 4  5  6
julia> b - U
ERROR: DimensionMismatch("matrix is not square: dimensions are (2, 3)")
Stacktrace:
[...]

Matrix factorizations

Matrix factorizations (a.k.a. matrix decompositions) compute the factorization of a matrix into a product of matrices, and are one of the central concepts in linear algebra.

The following table summarizes the types of matrix factorizations that have been implemented in Julia. Details of their associated methods can be found in the Standard Functions section of the Linear Algebra documentation.

Standard Functions

Linear algebra functions in Julia are largely implemented by calling functions from LAPACK. Sparse factorizations call functions from SuiteSparse.

Base.:* — Method.

*(A::AbstractMatrix, B::AbstractMatrix)

Matrix multiplication.

Examples

julia> [1 1; 0 1] * [1 0; 1 1]
2×2 Array{Int64,2}:
 2  1
 1  1

source

Base.:\ — Method.

\(A, B)

Matrix division using a polyalgorithm. For input matrices A and B, the result X is such that A*X == B when A is square. The solver that is used depends upon the structure of A. If A is upper or lower triangular (or diagonal), no factorization of A is required and the system is solved with either forward or backward substitution. For non-triangular square matrices, an LU factorization is used.

For rectangular A the result is the minimum-norm least squares solution computed by a pivoted QR factorization of A and a rank estimate of A 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.

Examples

julia> A = [1 0; 1 -2]; B = [32; -4];
julia> X = A \ B
2-element Array{Float64,1}:
 32.0
 18.0
julia> A * X == B
true

source

LinearAlgebra.dot — Function.

dot(x, y)
x ⋅ y

For any iterable containers x and y (including arrays of any dimension) of numbers (or any element type for which dot is defined), compute the dot product (or inner product or scalar product), i.e. the sum of dot(x[i],y[i]), as if they were vectors.

x ⋅ y (where ⋅ can be typed by tab-completing \cdot in the REPL) is a synonym for dot(x, y).

Examples

julia> dot(1:5, 2:6)
70
julia> x = fill(2., (5,5));
julia> y = fill(3., (5,5));
julia> dot(x, y)
150.0

source

dot(x, y)
x ⋅ y

Compute the dot product between two vectors. For complex vectors, the first vector is conjugated. When the vectors have equal lengths, calling dot is semantically equivalent to sum(dot(vx,vy) for (vx,vy) in zip(x, y)).

Examples

julia> dot([1; 1], [2; 3])
5
julia> dot([im; im], [1; 1])
0 - 2im

source

LinearAlgebra.cross — Function.

cross(x, y)
×(x,y)

Compute the cross product of two 3-vectors.

Examples

julia> a = [0;1;0]
3-element Array{Int64,1}:
 0
 1
 0
julia> b = [0;0;1]
3-element Array{Int64,1}:
 0
 0
 1
julia> cross(a,b)
3-element Array{Int64,1}:
 1
 0
 0

source

LinearAlgebra.factorize — Function.

factorize(A)

Compute a convenient factorization of A, based upon the type of the input matrix. factorize checks A to see if it is symmetric/triangular/etc. if A is passed as a generic matrix. factorize checks every element of A to verify/rule out each property. It will short-circuit as soon as it can rule out symmetry/triangular structure. The return value can be reused for efficient solving of multiple systems. For example: A=factorize(A); x=A\b; y=A\C.

If factorize is called on a Hermitian positive-definite matrix, for instance, then factorize will return a Cholesky factorization.

Examples

julia> A = Array(Bidiagonal(fill(1.0, (5, 5)), :U))
5×5 Array{Float64,2}:
 1.0  1.0  0.0  0.0  0.0
 0.0  1.0  1.0  0.0  0.0
 0.0  0.0  1.0  1.0  0.0
 0.0  0.0  0.0  1.0  1.0
 0.0  0.0  0.0  0.0  1.0
julia> factorize(A) # factorize will check to see that A is already factorized
5×5 Bidiagonal{Float64,Array{Float64,1}}:
 1.0  1.0   ⋅    ⋅    ⋅
  ⋅   1.0  1.0   ⋅    ⋅
  ⋅    ⋅   1.0  1.0   ⋅
  ⋅    ⋅    ⋅   1.0  1.0
  ⋅    ⋅    ⋅    ⋅   1.0

This returns a 5×5 Bidiagonal{Float64}, which can now be passed to other linear algebra functions (e.g. eigensolvers) which will use specialized methods for Bidiagonal types.

source

LinearAlgebra.Diagonal — Type.

Diagonal(A::AbstractMatrix)

Construct a matrix from the diagonal of A.

Examples

julia> A = [1 2 3; 4 5 6; 7 8 9]
3×3 Array{Int64,2}:
 1  2  3
 4  5  6
 7  8  9
julia> Diagonal(A)
3×3 Diagonal{Int64,Array{Int64,1}}:
 1  ⋅  ⋅
 ⋅  5  ⋅
 ⋅  ⋅  9

source

Diagonal(V::AbstractVector)

Construct a matrix with V as its diagonal.

Examples

julia> V = [1, 2]
2-element Array{Int64,1}:
 1
 2
julia> Diagonal(V)
2×2 Diagonal{Int64,Array{Int64,1}}:
 1  ⋅
 ⋅  2

source

LinearAlgebra.Bidiagonal — Type.

Bidiagonal(dv::V, ev::V, uplo::Symbol) where V <: AbstractVector

Constructs an upper (uplo=:U) or lower (uplo=:L) bidiagonal matrix using the given diagonal (dv) and off-diagonal (ev) vectors. The result is of type Bidiagonal and provides efficient specialized linear solvers, but may be converted into a regular matrix with convert(Array, _) (or Array(_) for short). The length of ev must be one less than the length of dv.

Examples

julia> dv = [1, 2, 3, 4]
4-element Array{Int64,1}:
 1
 2
 3
 4
julia> ev = [7, 8, 9]
3-element Array{Int64,1}:
 7
 8
 9
julia> Bu = Bidiagonal(dv, ev, :U) # ev is on the first superdiagonal
4×4 Bidiagonal{Int64,Array{Int64,1}}:
 1  7  ⋅  ⋅
 ⋅  2  8  ⋅
 ⋅  ⋅  3  9
 ⋅  ⋅  ⋅  4
julia> Bl = Bidiagonal(dv, ev, :L) # ev is on the first subdiagonal
4×4 Bidiagonal{Int64,Array{Int64,1}}:
 1  ⋅  ⋅  ⋅
 7  2  ⋅  ⋅
 ⋅  8  3  ⋅
 ⋅  ⋅  9  4

source

Bidiagonal(A, uplo::Symbol)

Construct a Bidiagonal matrix from the main diagonal of A and its first super- (if uplo=:U) or sub-diagonal (if uplo=:L).

Examples

julia> A = [1 1 1 1; 2 2 2 2; 3 3 3 3; 4 4 4 4]
4×4 Array{Int64,2}:
 1  1  1  1
 2  2  2  2
 3  3  3  3
 4  4  4  4
julia> Bidiagonal(A, :U) # contains the main diagonal and first superdiagonal of A
4×4 Bidiagonal{Int64,Array{Int64,1}}:
 1  1  ⋅  ⋅
 ⋅  2  2  ⋅
 ⋅  ⋅  3  3
 ⋅  ⋅  ⋅  4
julia> Bidiagonal(A, :L) # contains the main diagonal and first subdiagonal of A
4×4 Bidiagonal{Int64,Array{Int64,1}}:
 1  ⋅  ⋅  ⋅
 2  2  ⋅  ⋅
 ⋅  3  3  ⋅
 ⋅  ⋅  4  4

source

LinearAlgebra.SymTridiagonal — Type.

SymTridiagonal(dv::V, ev::V) where V <: AbstractVector

Construct a symmetric tridiagonal matrix from the diagonal (dv) and first sub/super-diagonal (ev), respectively. The result is of type SymTridiagonal and provides efficient specialized eigensolvers, but may be converted into a regular matrix with convert(Array, _) (or Array(_) for short).

Examples

julia> dv = [1, 2, 3, 4]
4-element Array{Int64,1}:
 1
 2
 3
 4
julia> ev = [7, 8, 9]
3-element Array{Int64,1}:
 7
 8
 9
julia> SymTridiagonal(dv, ev)
4×4 SymTridiagonal{Int64,Array{Int64,1}}:
 1  7  ⋅  ⋅
 7  2  8  ⋅
 ⋅  8  3  9
 ⋅  ⋅  9  4

source

SymTridiagonal(A::AbstractMatrix)

Construct a symmetric tridiagonal matrix from the diagonal and first sub/super-diagonal, of the symmetric matrix A.

Examples

julia> A = [1 2 3; 2 4 5; 3 5 6]
3×3 Array{Int64,2}:
 1  2  3
 2  4  5
 3  5  6
julia> SymTridiagonal(A)
3×3 SymTridiagonal{Int64,Array{Int64,1}}:
 1  2  ⋅
 2  4  5
 ⋅  5  6

source

LinearAlgebra.Tridiagonal — Type.

Tridiagonal(dl::V, d::V, du::V) where V <: AbstractVector

Construct a tridiagonal matrix from the first subdiagonal, diagonal, and first superdiagonal, respectively. The result is of type Tridiagonal and provides efficient specialized linear solvers, but may be converted into a regular matrix with convert(Array, _) (or Array(_) for short). The lengths of dl and du must be one less than the length of d.

Examples

julia> dl = [1, 2, 3];
julia> du = [4, 5, 6];
julia> d = [7, 8, 9, 0];
julia> Tridiagonal(dl, d, du)
4×4 Tridiagonal{Int64,Array{Int64,1}}:
 7  4  ⋅  ⋅
 1  8  5  ⋅
 ⋅  2  9  6
 ⋅  ⋅  3  0

source

Tridiagonal(A)

Construct a tridiagonal matrix from the first sub-diagonal, diagonal and first super-diagonal of the matrix A.

Examples

julia> A = [1 2 3 4; 1 2 3 4; 1 2 3 4; 1 2 3 4]
4×4 Array{Int64,2}:
 1  2  3  4
 1  2  3  4
 1  2  3  4
 1  2  3  4
julia> Tridiagonal(A)
4×4 Tridiagonal{Int64,Array{Int64,1}}:
 1  2  ⋅  ⋅
 1  2  3  ⋅
 ⋅  2  3  4
 ⋅  ⋅  3  4

source

LinearAlgebra.Symmetric — Type.

Symmetric(A, uplo=:U)

Construct a Symmetric view of the upper (if uplo = :U) or lower (if uplo = :L) triangle of the matrix A.

Examples

julia> A = [1 0 2 0 3; 0 4 0 5 0; 6 0 7 0 8; 0 9 0 1 0; 2 0 3 0 4]
5×5 Array{Int64,2}:
 1  0  2  0  3
 0  4  0  5  0
 6  0  7  0  8
 0  9  0  1  0
 2  0  3  0  4
julia> Supper = Symmetric(A)
5×5 Symmetric{Int64,Array{Int64,2}}:
 1  0  2  0  3
 0  4  0  5  0
 2  0  7  0  8
 0  5  0  1  0
 3  0  8  0  4
julia> Slower = Symmetric(A, :L)
5×5 Symmetric{Int64,Array{Int64,2}}:
 1  0  6  0  2
 0  4  0  9  0
 6  0  7  0  3
 0  9  0  1  0
 2  0  3  0  4

Note that Supper will not be equal to Slower unless A is itself symmetric (e.g. if A == transpose(A)).

source

LinearAlgebra.Hermitian — Type.

Hermitian(A, uplo=:U)

Construct a Hermitian view of the upper (if uplo = :U) or lower (if uplo = :L) triangle of the matrix A.

Examples

julia> A = [1 0 2+2im 0 3-3im; 0 4 0 5 0; 6-6im 0 7 0 8+8im; 0 9 0 1 0; 2+2im 0 3-3im 0 4];
julia> Hupper = Hermitian(A)
5×5 Hermitian{Complex{Int64},Array{Complex{Int64},2}}:
 1+0im  0+0im  2+2im  0+0im  3-3im
 0+0im  4+0im  0+0im  5+0im  0+0im
 2-2im  0+0im  7+0im  0+0im  8+8im
 0+0im  5+0im  0+0im  1+0im  0+0im
 3+3im  0+0im  8-8im  0+0im  4+0im
julia> Hlower = Hermitian(A, :L)
5×5 Hermitian{Complex{Int64},Array{Complex{Int64},2}}:
 1+0im  0+0im  6+6im  0+0im  2-2im
 0+0im  4+0im  0+0im  9+0im  0+0im
 6-6im  0+0im  7+0im  0+0im  3+3im
 0+0im  9+0im  0+0im  1+0im  0+0im
 2+2im  0+0im  3-3im  0+0im  4+0im

Note that Hupper will not be equal to Hlower unless A is itself Hermitian (e.g. if A == adjoint(A)).

All non-real parts of the diagonal will be ignored.

Hermitian(fill(complex(1,1), 1, 1)) == fill(1, 1, 1)

source

LinearAlgebra.LowerTriangular — Type.

LowerTriangular(A::AbstractMatrix)

Construct a LowerTriangular view of the matrix A.

Examples

julia> A = [1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0]
3×3 Array{Float64,2}:
 1.0  2.0  3.0
 4.0  5.0  6.0
 7.0  8.0  9.0
julia> LowerTriangular(A)
3×3 LowerTriangular{Float64,Array{Float64,2}}:
 1.0   ⋅    ⋅
 4.0  5.0   ⋅
 7.0  8.0  9.0

source

LinearAlgebra.UpperTriangular — Type.

UpperTriangular(A::AbstractMatrix)

Construct an UpperTriangular view of the matrix A.

Examples

julia> A = [1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0]
3×3 Array{Float64,2}:
 1.0  2.0  3.0
 4.0  5.0  6.0
 7.0  8.0  9.0
julia> UpperTriangular(A)
3×3 UpperTriangular{Float64,Array{Float64,2}}:
 1.0  2.0  3.0
  ⋅   5.0  6.0
  ⋅    ⋅   9.0

source

LinearAlgebra.UniformScaling — Type.

UniformScaling{T<:Number}

Generically sized uniform scaling operator defined as a scalar times the identity operator, λ*I. See also I.

Examples

julia> J = UniformScaling(2.)
UniformScaling{Float64}
2.0*I
julia> A = [1. 2.; 3. 4.]
2×2 Array{Float64,2}:
 1.0  2.0
 3.0  4.0
julia> J*A
2×2 Array{Float64,2}:
 2.0  4.0
 6.0  8.0

source

LinearAlgebra.lu — Function.

lu(A, pivot=Val(true); check = true) -> F::LU

Compute the LU factorization of A.

When check = true, an error is thrown if the decomposition fails. When check = false, responsibility for checking the decomposition’s validity (via issuccess) lies with the user.

In most cases, if A is a subtype S of AbstractMatrix{T} with an element type T supporting +, -, * and /, the return type is LU{T,S{T}}. If pivoting is chosen (default) the element type should also support abs and <.

The individual components of the factorization F can be accessed via getproperty:

Iterating the factorization produces the components F.L, F.U, and F.p.

The relationship between F and A is

F.L*F.U == A[F.p, :]

F further supports the following functions:

Examples

julia> A = [4 3; 6 3]
2×2 Array{Int64,2}:
 4  3
 6  3
julia> F = lu(A)
LU{Float64,Array{Float64,2}}
L factor:
2×2 Array{Float64,2}:
 1.0  0.0
 1.5  1.0
U factor:
2×2 Array{Float64,2}:
 4.0   3.0
 0.0  -1.5
julia> F.L * F.U == A[F.p, :]
true
julia> l, u, p = lu(A); # destructuring via iteration
julia> l == F.L && u == F.U && p == F.p
true

source

LinearAlgebra.lu! — Function.

lu!(A, pivot=Val(true); check = true) -> LU

lu! is the same as lu, but saves space by overwriting the input A, instead of creating a copy. An InexactError exception is thrown if the factorization produces a number not representable by the element type of A, e.g. for integer types.

Examples

julia> A = [4. 3.; 6. 3.]
2×2 Array{Float64,2}:
 4.0  3.0
 6.0  3.0
julia> F = lu!(A)
LU{Float64,Array{Float64,2}}
L factor:
2×2 Array{Float64,2}:
 1.0       0.0
 0.666667  1.0
U factor:
2×2 Array{Float64,2}:
 6.0  3.0
 0.0  1.0
julia> iA = [4 3; 6 3]
2×2 Array{Int64,2}:
 4  3
 6  3
julia> lu!(iA)
ERROR: InexactError: Int64(0.6666666666666666)
Stacktrace:
[...]

source

LinearAlgebra.cholesky — Function.

cholesky(A, Val(false); check = true) -> Cholesky

Compute the Cholesky factorization of a dense symmetric positive definite matrix A and return a Cholesky factorization. The matrix A can either be a Symmetric or Hermitian StridedMatrix or a perfectly symmetric or Hermitian StridedMatrix. The triangular Cholesky factor can be obtained from the factorization F with: F.L and F.U. The following functions are available for Cholesky objects: size, \, inv, det, logdet and isposdef.

When check = true, an error is thrown if the decomposition fails. When check = false, responsibility for checking the decomposition’s validity (via issuccess) lies with the user.

Examples

julia> A = [4. 12. -16.; 12. 37. -43.; -16. -43. 98.]
3×3 Array{Float64,2}:
   4.0   12.0  -16.0
  12.0   37.0  -43.0
 -16.0  -43.0   98.0
julia> C = cholesky(A)
Cholesky{Float64,Array{Float64,2}}
U factor:
3×3 UpperTriangular{Float64,Array{Float64,2}}:
 2.0  6.0  -8.0
  ⋅   1.0   5.0
  ⋅    ⋅    3.0
julia> C.U
3×3 UpperTriangular{Float64,Array{Float64,2}}:
 2.0  6.0  -8.0
  ⋅   1.0   5.0
  ⋅    ⋅    3.0
julia> C.L
3×3 LowerTriangular{Float64,Array{Float64,2}}:
  2.0   ⋅    ⋅
  6.0  1.0   ⋅
 -8.0  5.0  3.0
julia> C.L * C.U == A
true

source

cholesky(A, Val(true); tol = 0.0, check = true) -> CholeskyPivoted

Compute the pivoted Cholesky factorization of a dense symmetric positive semi-definite matrix A and return a CholeskyPivoted factorization. The matrix A can either be a Symmetric or Hermitian StridedMatrix or a perfectly symmetric or Hermitian StridedMatrix. The triangular Cholesky factor can be obtained from the factorization F with: F.L and F.U. The following functions are available for CholeskyPivoted objects: size, \, inv, det, and rank. The argument tol determines the tolerance for determining the rank. For negative values, the tolerance is the machine precision.

When check = true, an error is thrown if the decomposition fails. When check = false, responsibility for checking the decomposition’s validity (via issuccess) lies with the user.

source

LinearAlgebra.cholesky! — Function.

cholesky!(A, Val(false); check = true) -> Cholesky

The same as cholesky, but saves space by overwriting the input A, instead of creating a copy. An InexactError exception is thrown if the factorization produces a number not representable by the element type of A, e.g. for integer types.

Examples

julia> A = [1 2; 2 50]
2×2 Array{Int64,2}:
 1   2
 2  50
julia> cholesky!(A)
ERROR: InexactError: Int64(6.782329983125268)
Stacktrace:
[...]

source

cholesky!(A, Val(true); tol = 0.0, check = true) -> CholeskyPivoted

The same as cholesky, but saves space by overwriting the input A, instead of creating a copy. An InexactError exception is thrown if the factorization produces a number not representable by the element type of A, e.g. for integer types.

source

LinearAlgebra.lowrankupdate — Function.

lowrankupdate(C::Cholesky, v::StridedVector) -> CC::Cholesky

Update a Cholesky factorization C with the vector v. If A = C.U'C.U then CC = cholesky(C.U'C.U + v*v') but the computation of CC only uses O(n^2) operations.

source

LinearAlgebra.lowrankdowndate — Function.

lowrankdowndate(C::Cholesky, v::StridedVector) -> CC::Cholesky

Downdate a Cholesky factorization C with the vector v. If A = C.U'C.U then CC = cholesky(C.U'C.U - v*v') but the computation of CC only uses O(n^2) operations.

source

LinearAlgebra.lowrankupdate! — Function.

lowrankupdate!(C::Cholesky, v::StridedVector) -> CC::Cholesky

Update a Cholesky factorization C with the vector v. If A = C.U'C.U then CC = cholesky(C.U'C.U + v*v') but the computation of CC only uses O(n^2) operations. The input factorization C is updated in place such that on exit C == CC. The vector v is destroyed during the computation.

source

LinearAlgebra.lowrankdowndate! — Function.

lowrankdowndate!(C::Cholesky, v::StridedVector) -> CC::Cholesky

Downdate a Cholesky factorization C with the vector v. If A = C.U'C.U then CC = cholesky(C.U'C.U - v*v') but the computation of CC only uses O(n^2) operations. The input factorization C is updated in place such that on exit C == CC. The vector v is destroyed during the computation.

source

LinearAlgebra.ldlt — Function.

ldlt(S::SymTridiagonal) -> LDLt

Compute an LDLt factorization of the real symmetric tridiagonal matrix S such that S = L*Diagonal(d)*L' where L is a unit lower triangular matrix and d is a vector. The main use of an LDLt factorization F = ldlt(S) is to solve the linear system of equations Sx = b with F\b.

Examples

julia> S = SymTridiagonal([3., 4., 5.], [1., 2.])
3×3 SymTridiagonal{Float64,Array{Float64,1}}:
 3.0  1.0   ⋅
 1.0  4.0  2.0
  ⋅   2.0  5.0
julia> ldltS = ldlt(S);
julia> b = [6., 7., 8.];
julia> ldltS \ b
3-element Array{Float64,1}:
 1.7906976744186047
 0.627906976744186
 1.3488372093023255
julia> S \ b
3-element Array{Float64,1}:
 1.7906976744186047
 0.627906976744186
 1.3488372093023255

source

LinearAlgebra.ldlt! — Function.

ldlt!(S::SymTridiagonal) -> LDLt

Same as ldlt, but saves space by overwriting the input S, instead of creating a copy.

Examples

julia> S = SymTridiagonal([3., 4., 5.], [1., 2.])
3×3 SymTridiagonal{Float64,Array{Float64,1}}:
 3.0  1.0   ⋅
 1.0  4.0  2.0
  ⋅   2.0  5.0
julia> ldltS = ldlt!(S);
julia> ldltS === S
false
julia> S
3×3 SymTridiagonal{Float64,Array{Float64,1}}:
 3.0       0.333333   ⋅
 0.333333  3.66667   0.545455
  ⋅        0.545455  3.90909

source

LinearAlgebra.qr — Function.

qr(A, pivot=Val(false)) -> F

Compute the QR factorization of the matrix A: an orthogonal (or unitary if A is complex-valued) matrix Q, and an upper triangular matrix R such that

\[A = Q R\]

The returned object F stores the factorization in a packed format:

• if pivot == Val(true) then F is a QRPivoted object,

• otherwise if the element type of A is a BLAS type (Float32, Float64, ComplexF32 or ComplexF64), then F is a QRCompactWY object,

• otherwise F is a QR object.

The individual components of the decomposition F can be retrieved via property accessors:

• F.Q: the orthogonal/unitary matrix Q
• F.R: the upper triangular matrix R
• F.p: the permutation vector of the pivot (QRPivoted only)
• F.P: the permutation matrix of the pivot (QRPivoted only)

Iterating the decomposition produces the components Q, R, and if extant p.

The following functions are available for the QR objects: inv, size, and \. When A is rectangular, \ will return a least squares solution and if the solution is not unique, the one with smallest norm is returned. When A 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/square Q is allowed, i.e. both F.Q*F.R and F.Q*A are supported. A Q matrix can be converted into a regular matrix with Matrix. This operation returns the “thin” Q factor, i.e., if A is m×n with m>=n, then Matrix(F.Q) yields an m×n matrix with orthonormal columns. To retrieve the “full” Q factor, an m×m orthogonal matrix, use F.Q*Matrix(I,m,m). If m<=n, then Matrix(F.Q) yields an m×m orthogonal matrix.

Examples

julia> A = [3.0 -6.0; 4.0 -8.0; 0.0 1.0]
3×2 Array{Float64,2}:
 3.0  -6.0
 4.0  -8.0
 0.0   1.0
julia> F = qr(A)
LinearAlgebra.QRCompactWY{Float64,Array{Float64,2}}
Q factor:
3×3 LinearAlgebra.QRCompactWYQ{Float64,Array{Float64,2}}:
 -0.6   0.0   0.8
 -0.8   0.0  -0.6
  0.0  -1.0   0.0
R factor:
2×2 Array{Float64,2}:
 -5.0  10.0
  0.0  -1.0
julia> F.Q * F.R == A
true

Note

qr returns multiple types because LAPACK uses several representations that minimize the memory storage requirements of products of Householder elementary reflectors, so that the Q and R matrices can be stored compactly rather as two separate dense matrices.

source

LinearAlgebra.qr! — Function.

qr!(A, pivot=Val(false))

qr! is the same as qr when A is a subtype of StridedMatrix, but saves space by overwriting the input A, instead of creating a copy. An InexactError exception is thrown if the factorization produces a number not representable by the element type of A, e.g. for integer types.

Examples

julia> a = [1. 2.; 3. 4.]
2×2 Array{Float64,2}:
 1.0  2.0
 3.0  4.0
julia> qr!(a)
LinearAlgebra.QRCompactWY{Float64,Array{Float64,2}}
Q factor:
2×2 LinearAlgebra.QRCompactWYQ{Float64,Array{Float64,2}}:
 -0.316228  -0.948683
 -0.948683   0.316228
R factor:
2×2 Array{Float64,2}:
 -3.16228  -4.42719
  0.0      -0.632456
julia> a = [1 2; 3 4]
2×2 Array{Int64,2}:
 1  2
 3  4
julia> qr!(a)
ERROR: InexactError: Int64(-3.1622776601683795)
Stacktrace:
[...]

source

LinearAlgebra.QR — Type.

QR <: Factorization

A QR matrix factorization stored in a packed format, typically obtained from qr. If $A$ is an m×n matrix, then

\[A = Q R\]

where $Q$ is an orthogonal/unitary matrix and $R$ is upper triangular. The matrix $Q$ is stored as a sequence of Householder reflectors $v_i$ and coefficients $\tau_i$ where:

\[Q = \prod_{i=1}^{\min(m,n)} (I - \tau_i v_i v_i^T).\]

Iterating the decomposition produces the components Q and R.

The object has two fields:

• factors is an m×n matrix.

  • The upper triangular part contains the elements of $R$, that is R = triu(F.factors) for a QR object F.

  • The subdiagonal part contains the reflectors $v_i$ stored in a packed format where $v_i$ is the $i$th column of the matrix V = I + tril(F.factors, -1).

• τ is a vector of length min(m,n) containing the coefficients $au_i$.

source

LinearAlgebra.QRCompactWY — Type.

QRCompactWY <: Factorization

A QR matrix factorization stored in a compact blocked format, typically obtained from qr. If $A$ is an m×n matrix, then

\[A = Q R\]

where $Q$ is an orthogonal/unitary matrix and $R$ is upper triangular. It is similar to the QR format except that the orthogonal/unitary matrix $Q$ is stored in Compact WY format [Schreiber1989], as a lower trapezoidal matrix $V$ and an upper triangular matrix $T$ where

\[Q = \prod_{i=1}^{\min(m,n)} (I - \tau_i v_i v_i^T) = I - V T V^T\]

such that $v_i$ is the $i$th column of $V$, and $au_i$ is the $i$th diagonal element of $T$.

Iterating the decomposition produces the components Q and R.

The object has two fields:

• factors, as in the QR type, is an m×n matrix.

  • The upper triangular part contains the elements of $R$, that is R = triu(F.factors) for a QR object F.

  • The subdiagonal part contains the reflectors $v_i$ stored in a packed format such that V = I + tril(F.factors, -1).

• T is a square matrix with min(m,n) columns, whose upper triangular part gives the matrix $T$ above (the subdiagonal elements are ignored).

Note

This format should not to be confused with the older WY representation [Bischof1987].

[Bischof1987]

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

source

LinearAlgebra.QRPivoted — Type.

QRPivoted <: Factorization

A QR matrix factorization with column pivoting in a packed format, typically obtained from qr. If $A$ is an m×n matrix, then

\[A P = Q R\]

where $P$ is a permutation matrix, $Q$ is an orthogonal/unitary matrix and $R$ is upper triangular. The matrix $Q$ is stored as a sequence of Householder reflectors:

\[Q = \prod_{i=1}^{\min(m,n)} (I - \tau_i v_i v_i^T).\]

Iterating the decomposition produces the components Q, R, and p.

The object has three fields:

• factors is an m×n matrix.

  • The upper triangular part contains the elements of $R$, that is R = triu(F.factors) for a QR object F.

  • The subdiagonal part contains the reflectors $v_i$ stored in a packed format where $v_i$ is the $i$th column of the matrix V = I + tril(F.factors, -1).

• τ is a vector of length min(m,n) containing the coefficients $au_i$.

• jpvt is an integer vector of length n corresponding to the permutation $P$.

source

LinearAlgebra.lq! — Function.

lq!(A) -> LQ

Compute the LQ factorization of A, using the input matrix as a workspace. See also lq.

source

LinearAlgebra.lq — Function.

lq(A) -> S::LQ

Compute the LQ decomposition of A. The decomposition’s lower triangular component can be obtained from the LQ object S via S.L, and the orthogonal/unitary component via S.Q, such that A ≈ S.L*S.Q.

Iterating the decomposition produces the components S.L and S.Q.

The LQ decomposition is the QR decomposition of transpose(A).

Examples

julia> A = [5. 7.; -2. -4.]
2×2 Array{Float64,2}:
  5.0   7.0
 -2.0  -4.0
julia> S = lq(A)
LQ{Float64,Array{Float64,2}} with factors L and Q:
[-8.60233 0.0; 4.41741 -0.697486]
[-0.581238 -0.813733; -0.813733 0.581238]
julia> S.L * S.Q
2×2 Array{Float64,2}:
  5.0   7.0
 -2.0  -4.0
julia> l, q = S; # destructuring via iteration
julia> l == S.L &&  q == S.Q
true

source

LinearAlgebra.bunchkaufman — Function.

bunchkaufman(A, rook::Bool=false; check = true) -> S::BunchKaufman

Compute the Bunch-Kaufman [Bunch1977] factorization of a Symmetric or Hermitian matrix A as $P’*U*D*U’*P$ or $P’*L*D*L’*P$, depending on which triangle is stored in A, and return a BunchKaufman object. Note that if A is complex symmetric then U' and L' denote the unconjugated transposes, i.e. transpose(U) and transpose(L).

Iterating the decomposition produces the components S.D, S.U or S.L as appropriate given S.uplo, and S.p.

If rook is true, rook pivoting is used. If rook is false, rook pivoting is not used.

When check = true, an error is thrown if the decomposition fails. When check = false, responsibility for checking the decomposition’s validity (via issuccess) lies with the user.

The following functions are available for BunchKaufman objects: size, \, inv, issymmetric, ishermitian, getindex.

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

Examples

julia> A = [1 2; 2 3]
2×2 Array{Int64,2}:
 1  2
 2  3
julia> S = bunchkaufman(A)
BunchKaufman{Float64,Array{Float64,2}}
D factor:
2×2 Tridiagonal{Float64,Array{Float64,1}}:
 -0.333333  0.0
  0.0       3.0
U factor:
2×2 UnitUpperTriangular{Float64,Array{Float64,2}}:
 1.0  0.666667
  ⋅   1.0
permutation:
2-element Array{Int64,1}:
 1
 2
julia> d, u, p = S; # destructuring via iteration
julia> d == S.D && u == S.U && p == S.p
true

source

LinearAlgebra.bunchkaufman! — Function.

bunchkaufman!(A, rook::Bool=false; check = true) -> BunchKaufman

bunchkaufman! is the same as bunchkaufman, but saves space by overwriting the input A, instead of creating a copy.

source

LinearAlgebra.eigvals — Function.

eigvals(A; permute::Bool=true, scale::Bool=true) -> values

Return the eigenvalues of A.

For general non-symmetric matrices it is possible to specify how the matrix is balanced before the eigenvalue calculation. The option permute=true permutes the matrix to become closer to upper triangular, and scale=true scales the matrix by its diagonal elements to make rows and columns more equal in norm. The default is true for both options.

Examples

julia> diag_matrix = [1 0; 0 4]
2×2 Array{Int64,2}:
 1  0
 0  4
julia> eigvals(diag_matrix)
2-element Array{Float64,1}:
 1.0
 4.0