decomposition

Matrix decomposition for solving linear systems

Description

decomposition creates reusable matrix decompositions (LU, LDL, Cholesky, QR, and more) that enable you to solve linear systems (Ax = b or xA = b) more efficiently. For example, after computing dA = decomposition(A) the call dA\b returns the same vector as A\b, but is typically much faster. decomposition objects are well-suited to solving problems that require repeated solutions, since the decomposition of the coefficient matrix does not need to be performed multiple times.

You can use a decomposition object dA with many of the same operators you might use on the original coefficient matrix A:

• Complex conjugate transpose dA'

• Negation -dA

• Multiply or divide by a scalar using c*dA or dA/c.

• Solve a linear system Ax = b using x = dA\b.

• Solve a linear system xA = b using x = b/dA.

Creation

Syntax

dA = decomposition(A)

dA = decomposition(A,type)

dA = decomposition(A,type,triangularFlag)

dA = decomposition(___,Name,Value)

Description

dA = decomposition(A) returns a decomposition of matrix A that you can use to solve linear systems more efficiently. The decomposition type is automatically chosen based on the properties of the input matrix.

dA = decomposition(A,type) specifies the type of decomposition to perform. type can be 'qr', 'cod', 'lu', 'ldl', 'chol', 'triangular', 'permutedTriangular', 'banded', 'hessenberg', or 'diagonal'.

dA = decomposition(A,type,triangularFlag) specifies that only the upper or lower triangular portion of A is to be used in the decomposition. triangularFlag can be 'upper' or 'lower'. With this syntax, the decomposition type must be 'ldl', 'chol', or 'triangular'.

dA = decomposition(___,Name,Value) specifies additional options using one or more Name,Value pair arguments using any of the previous syntaxes. For example, dA = decomposition(A,'CheckCondition',false) specifies not to throw a warning based on the condition of A while solving dA\b.

Input Arguments

A — Coefficient matrix

matrix

Coefficient matrix. The coefficient matrix appears in the system of linear equations on the left as Ax = b or on the right as xA = b.

Data Types: single | double
Complex Number Support: Yes

type — Decomposition type

'auto' (default) | 'qr' | 'cod' | 'lu' | 'ldl' | 'chol' | 'triangular' | 'permutedTriangular' | 'banded' | 'hessenberg' | 'diagonal'

Decomposition type, specified as one of the options in these tables.

These options work for any coefficient matrix.

Value	Matrix Decomposition of A	Notes
'auto' (default)	N/A	Automatic selection of the matrix decomposition type based on the properties of the coefficient matrix. For information about how the decomposition is chosen, see the Algorithms section of mldivide.
'qr'	$$AP=QR$$ Q is unitary, R is upper triangular, and P is a permutation matrix.	QR decompositions give a least-squares solution. If type is 'qr', then you cannot solve A'\B or B/A. Instead, use 'cod' for problems with those forms.
'cod'	$$A=QR{Z}^{*}$$ R is upper triangular, and both Q and Z have orthonormal columns.	Complete orthogonal decompositions give a minimum norm least-squares solution.

For square coefficient matrices, you also can use these options.

Value	Matrix Decomposition of A	Notes
'lu'	Dense matrices: $$PA=LU$$ L is lower triangular, U is upper triangular, and P is a permutation matrix. Sparse matrices: $$P\left(R\backslash A\right)Q=LU$$ P and Q are permutation matrices and R is a diagonal scaling matrix.
'ldl'	Dense matrices: $$P^{*}AP=LD{L}^{*}$$ L is a lower triangular matrix with 1s on the diagonal, D is a diagonal matrix, and P is a permutation matrix. Sparse matrices: $$P^{*}SASP=LD{L}^{*}$$ S is a scaling matrix.	A must be symmetric.
'chol'	Dense matrices: $$A=L{L}^{*}$$ L is lower triangular. Sparse matrices: $$A=PL{L}^{*}{P}^{*}$$ P is a permutation matrix.	A must be symmetric positive definite.
'triangular'	$$A=T$$ T is triangular.	A must be triangular.
'permutedTriangular'	$$A=PT{Q}^{*}$$ T is triangular, and both P and Q are permutation matrices.	A must be a permutation of a triangular matrix.
'banded'	$$A=P^{*}LU$$ P is a permutation matrix and both L and U are banded.	Most effective for matrices with a low bandwidth. See bandwidth for more information.
'hessenberg'	$$A=P^{*}LU$$ P is a permutation matrix, L is banded, and U is upper triangular.	A must have zeros below the first subdiagonal.
'diagonal'	$$A=D$$ D is diagonal.	A must be diagonal.

triangularFlag — Flag to use only upper or lower triangular portion of coefficient matrix

'upper' | 'lower'

Flag to use only upper or lower triangular portion of coefficient matrix, specified as either 'upper' or 'lower'. This option supports the 'triangular', 'chol', and 'ldl' decomposition types.

• 'triangular' — If both an upper and lower triangular matrix are stored in the same matrix, then use triangularFlag to specify only one of the triangular matrices.

• 'chol' and 'ldl' — Use triangularFlag to avoid symmetrizing a nearly symmetric coefficient matrix.

Name-Value Arguments

Specify optional pairs of arguments as Name1=Value1,...,NameN=ValueN, where Name is the argument name and Value is the corresponding value. Name-value arguments must appear after other arguments, but the order of the pairs does not matter.

Before R2021a, use commas to separate each name and value, and enclose Name in quotes.

Example: dA = decomposition(A,'qr','CheckCondition',false) performs a QR decomposition of A and turns off warnings about the condition of the coefficient matrix when it is used to solve a linear system.

General Parameters

CheckCondition — Toggle to check condition of coefficient matrix

true (default) | false

Toggle to check condition of coefficient matrix, specified as the comma-separated pair consisting of 'CheckCondition' and either logical 1 (true) or logical 0 (false). If CheckCondition is true and the coefficient matrix is badly conditioned or of low rank, then solving linear systems using mldivide (\) or mrdivide (/) produces warnings.

Data Types: logical

RankTolerance — Rank tolerance

nonnegative scalar

Rank tolerance, specified as a nonnegative scalar. Specifying the tolerance can help prevent the solution from being susceptible to random noise in the coefficient matrix.

decomposition computes the rank of A as the number of diagonal elements in the R matrix of the QR decomposition [Q,R,p] = qr(A,0) with absolute value larger than tol. If the rank of A is k, then a low-rank approximation of A is formed by multiplying the first k columns of Q by the first k rows of R. Changing the tolerance affects this low-rank approximation of A.

Note

This option only applies when 'Type' is 'qr' or 'cod', or when 'Type' is 'auto' and A is rectangular. Otherwise, this option is ignored.

Sparse Matrix Parameters

BandDensity — Band density threshold

0.5 (default) | scalar

Band density threshold, specified as a scalar value in the range [0 1]. The value of 'BandDensity' determines how dense a sparse, banded coefficient matrix must be for the banded solver to be used by mldivide (\) or mrdivide (/) when solving a system of equations. If the band density of the coefficient matrix is larger than the specified band density, then the banded solver is used.

The band density is defined as: (# nonzeros in the band) / (# elements in the band). A value of 1.0 indicates to never use the banded solver.

LDLPivotTolerance — Pivot tolerance for LDL factorization

0.01 (default) | scalar

Pivot tolerance for LDL factorization, specified as a scalar value in the interval [0 0.5]. Using smaller values of the pivot tolerance can give faster factorization times and fewer entries, but also can result in a less stable factorization.

This pivot tolerance is the same that ldl uses for real sparse matrices.

LUPivotTolerance — Pivot tolerance for LU factorization

[0.1 0.001] (default) | scalar | vector

Pivot tolerance for LU factorization, specified as a scalar or vector. Specify a scalar value to change the first element in the tolerance vector, or specify a two-element vector to change both values. Smaller pivot tolerances tend to lead to sparser LU factors, but the solution can become inaccurate. Larger values can lead to a more accurate solution, but not always, and usually increase the total work and memory usage.

This pivot tolerance is the same that lu uses for sparse matrices.

Properties

MatrixSize — Size of coefficient matrix

vector

This property is read-only.

Size of coefficient matrix, returned as a two-element row vector.

Data Types: double

Type — Decomposition type

'qr' | 'cod' | 'lu' | 'ldl' | 'chol' | 'triangular' | 'permutedTriangular' | 'banded' | 'hessenberg' | 'diagonal'

This property is read-only.

Decomposition type, returned as 'qr', 'cod', 'lu', 'ldl', 'chol', 'triangular', 'permutedTriangular', 'banded', 'hessenberg', or 'diagonal'.

Data Types: char

CheckCondition — Toggle to check condition of coefficient matrix

true (default) | false

Toggle to check condition of coefficient matrix, specified as either logical 1 (true) or logical 0 (false). If CheckCondition is true and the coefficient matrix is badly conditioned or of low rank, then solving linear systems using mldivide (\) or mrdivide (/) produces warnings.

Data Types: logical

Datatype — Data type of coefficient matrix

'double' | 'single'

This property is read-only.

Data type of coefficient matrix, returned as either 'double' or 'single'.

Data Types: char

IsConjugateTransposed — Indicator that coefficient matrix is complex conjugate transposed

false (default) | true

This property is read-only.

Indicator that coefficient matrix is complex conjugate transposed, returned as either logical 1 (true) or logical 0 (false). This indicator is false by default for any decomposition object constructed from the coefficient matrix. However, the value is true if you use the ctranspose operator on a decomposition object in an expression, such as dA'\b. In that case, dA' is the same decomposition object as dA, but with a value of true for IsConjugateTransposed.

Data Types: logical

IsReal — Indicator that coefficient matrix is real

true | false

This property is read-only.

Indicator that coefficient matrix is real, returned as either logical 1 (true) or logical 0 (false). A value of false indicates that the coefficient matrix contains complex numbers.

Data Types: logical

IsSparse — Indicator that coefficient matrix is sparse

true | false

This property is read-only.

Indicator that coefficient matrix is sparse, returned as either logical 1 (true) or logical 0 (false).

Data Types: logical

ScaleFactor — Multiplicative scale factor for coefficient matrix

1 (default) | scalar

This property is read-only.

Multiplicative scale factor for coefficient matrix, returned as a scalar. The default value of 1 indicates that the coefficient matrix is not scaled. However, when you multiply or divide the decomposition object by a scalar, the value of ScaleFactor changes. For example, 3*dA is a decomposition object equivalent to dA, but with a value of 3 for ScaleFactor.

Data Types: double
Complex Number Support: Yes

Object Functions

The primary functions and operators that you can use with decomposition objects are related to solving linear systems of equations. If the decomposition type is 'qr', then you cannot solve A'\B or B/A. Instead, use 'cod' for problems with those forms.

ctranspose	Complex conjugate transpose
mldivide	Solve systems of linear equations Ax = B for x
mrdivide	Solve systems of linear equations xA = B for x
isIllConditioned	Determine whether matrix is ill conditioned

You also can check the condition number or rank of the underlying matrix of decomposition objects. Since different algorithms are employed, the results of using these functions on the decomposition object can differ compared to using the same functions directly on the coefficient matrix.

rank	Only the 1-input form rank(dA) is supported. The decomposition type must be 'qr' or 'cod'. The value of the rank depends on the choice of RankTolerance, if specified.
rcond	Runs the same condition check that backslash \ uses to determine whether to issue a warning. Supports all decomposition types except for 'qr' and 'cod'.

Examples

Solve Linear System with Several Right-Hand Sides

Show how using decomposition objects can improve the efficiency of solving $\mathrm{Ax}=\mathit{b}$ with many right-hand sides.

The inverse iteration is an iterative eigenvalue algorithm that solves linear systems with many right-hand sides. It is a method to iteratively compute an eigenvalue of a matrix starting from a guess of the corresponding eigenvector. Each iteration computes x = A\x, and then scales x by its norm.

Create a sparse matrix A and random starting vectors x1 and x2.

n = 1e3;
rng default % for reproducibility
A = sprandn(n,n,0.2) + speye(n);
x1 = randn(n,1);
x2 = x1;

Apply 100 iterations of the inverse iteration algorithm using backslash to calculate an eigenvalue of A.

tic
for ii=1:100
    x1 = A \ x1;
    x1 = x1 / norm(x1);
end
toc

Elapsed time is 11.755453 seconds.

lambda = x1'*A*x1

lambda = 
-0.6707

Now use a decomposition object to solve the same problem.

tic
dA = decomposition(A); 
for ii=1:100
    x2 = dA \ x2;
    x2 = x2 / norm(x2);
end
toc

Elapsed time is 1.327210 seconds.

lambda = x2'*A*x2

lambda = 
-0.6707

The performance of the algorithm improves dramatically because the matrix A does not need to be factorized during each iteration. Also, even though the backslash algorithm can be improved by performing an LU decomposition of A before the for-loop, the decomposition object gives access to all of the same performance gains without requiring that you write complex code.

Select Decomposition Type

Choose a decomposition type to override the automatic default selection based on the input matrix.

Create a coefficient matrix and decompose the matrix using the default selection of decomposition type.

A = ones(3);
dA = decomposition(A)

dA = 
  decomposition with properties:

    MatrixSize: [3 3]
          Type: 'lu'

  Show all properties

Solve the linear system using a vector of ones for the right-hand side.

b = ones(3,1);
x = dA\b

Warning: Matrix is singular to working precision.

x = 3×1

   NaN
   NaN
   NaN

Specify the decomposition type to use the 'qr' method instead of the default 'ldl' method. This forces backslash (\) to find a least-squares solution to the problem instead of returning a vector of NaNs.

dA_qr = decomposition(A,'qr')

dA_qr = 
  decomposition with properties:

    MatrixSize: [3 3]
          Type: 'qr'

  Show all properties

x = dA_qr\b

Warning: Rank deficient, rank = 1, tol =  1.153778e-15.

x = 3×1

     1
     0
     0

Use Triangular Portion of Matrix

Specify 'upper' to use only the upper triangular portion of an input matrix in the decomposition.

Create a coefficient matrix. Construct a triangular decomposition for the matrix using only the upper triangular portion. This option can be useful in cases where both an upper triangular and lower triangular matrix are stored in the same matrix.

A = randi([0 5],10)

A = 10×10

     4     0     3     4     2     1     4     5     2     0
     5     5     0     0     2     4     1     1     4     0
     0     5     5     1     4     3     3     4     3     3
     5     2     5     0     4     0     4     1     3     4
     3     4     4     0     1     0     5     5     5     5
     0     0     4     4     2     2     5     2     1     0
     1     2     4     4     2     5     3     1     4     3
     3     5     2     1     3     2     0     1     4     2
     5     4     3     5     4     3     0     3     2     0
     5     5     1     0     4     1     1     2     3     2

dA = decomposition(A,'triangular','upper')

dA = 
  decomposition with properties:

    MatrixSize: [10 10]
          Type: 'triangular'

  Show all properties

Turn Off Matrix Condition Warnings

Use the 'CheckCondition' name-value pair to turn off warnings based on the condition of the coefficient matrix when solving a linear system using decomposition.

Create a coefficient matrix that is ill conditioned. In this matrix, averaging together the first two columns produces the third column.

A = [1 2 1.5; 3 4 3.5; 5 6 5.5]

A = 3×3

    1.0000    2.0000    1.5000
    3.0000    4.0000    3.5000
    5.0000    6.0000    5.5000

Solve a linear system $\mathrm{Ax}=\mathit{b}$ using a vector of 1s for the right-hand side. mldivide produces a warning about the conditioning of the coefficient matrix.

b = ones(3,1);
x = A\b

Warning: Matrix is singular to working precision.

x = 3×1

   NaN
   Inf
  -Inf

Now create a decomposition object for the matrix and solve the same linear system. Specify 'CheckCondition' as false so that mldivide does not check the condition of the coefficient matrix. Even though the same solution is returned, mldivide does not display the warning message.

dA = decomposition(A,'CheckCondition',false);
x = dA\b

x = 3×1

   NaN
   Inf
  -Inf

Use the isIllConditioned function to check whether the decomposition object is based on an ill-conditioned matrix.

tf = isIllConditioned(dA)

tf = logical
   1

References

[1] Davis, Timothy A. “Algorithm 930: FACTORIZE: An Object-Oriented Linear System Solver for MATLAB.” ACM Transactions on Mathematical Software 39, no. 4 (July 2013): 1–18. https://doi.org/10.1145/2491491.2491498.

Extended Capabilities

Thread-Based Environment

Run code in the background using MATLAB® backgroundPool or accelerate code with Parallel Computing Toolbox™ ThreadPool.

This function fully supports thread-based environments. For more information, see Run MATLAB Functions in Thread-Based Environment.

Distributed Arrays

Partition large arrays across the combined memory of your cluster using Parallel Computing Toolbox™.

Usage notes and limitations:

• If A is dense, then the values 'ldl', 'cod', and 'hessenberg' are not supported for type.

• If A is sparse, then the values 'chol', 'cod', and 'hessenberg' are not supported for type.

• If A is sparse and type is 'lu', then CheckCondition must be false, as rcond is unable to determine the reciprocal condition number.

• If A is sparse, rcond(decomposition(A,'lu','CheckCondition',false)), MATLAB^® returns [].

For more information, see Run MATLAB Functions with Distributed Arrays (Parallel Computing Toolbox).

Version History

Introduced in R2017b