decomposition - Matrix decomposition for solving linear systems - MATLAB (original) (raw)

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 matrixA:

• Complex conjugate transpose dA'
• Negation -dA
• Multiply or divide by a scalar using c*dA ordA/c.
• Solve a linear system Ax = b using x = dA\b.
• Solve a linear system xA = b using x = b/dA.

Creation

Syntax

Description

`dA` = decomposition([A](#d126e387352)) 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.

example

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

example

`dA` = decomposition([A](#d126e387352),[type](#d126e387387),[triangularFlag](#d126e387755)) specifies that only the upper or lower triangular portion ofA 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'.

example

`dA` = decomposition(___,[Name,Value](#namevaluepairarguments)) 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 solvingdA\b.

example

Input Arguments

expand all

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.

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

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 asName1=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

expand all

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 logical1 (true) or logical0 (false). IfCheckCondition is true and the coefficient matrix is badly conditioned or of low rank, then solving linear systems using mldivide (\) ormrdivide (/) 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 ofA as the number of diagonal elements in theR matrix of the QR decomposition[Q,R,p] = qr(A,0) with absolute value larger than tol. If the rank of A isk, then a low-rank approximation ofA is formed by multiplying the firstk columns of Q by the first k rows of R. Changing the tolerance affects this low-rank approximation ofA.

Note

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

Sparse Matrix Parameters

expand all

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

expand all

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 logical0 (false). IfCheckCondition 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 isfalse by default for anydecomposition object constructed from the coefficient matrix. However, the value is true if you use thectranspose operator on adecomposition object in an expression, such asdA'\b. In that case, dA' is the same decomposition object as dA, but with a value of true forIsConjugateTransposed.

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 logical1 (true) or logical0 (false). A value offalse 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 logical1 (true) or logical0 (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 thedecomposition object by a scalar, the value ofScaleFactor changes. For example,3*dA is a decomposition object equivalent to dA, but with a value of3 for ScaleFactor.

Data Types: double
Complex Number Support: Yes

Object Functions

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

You also can check the condition number or rank of the underlying matrix ofdecomposition 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.

Examples

collapse all

Solve Linear System with Several Right-Hand Sides

Show how using decomposition objects can improve the efficiency of solving Ax=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 12.394779 seconds.

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.321645 seconds.

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.

Warning: Matrix is singular to working precision.

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

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

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 = 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 Ax=b using a vector of 1s for the right-hand side. mldivide produces a warning about the conditioning of the coefficient matrix.

Warning: Matrix is singular to working precision.

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

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

tf = isIllConditioned(dA)

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 fortype.
• If A is sparse, then the values'chol', 'cod', and'hessenberg' are not supported fortype.
• If A is sparse and type is 'lu', then CheckCondition must be false, asrcond 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