Code covered by the BSD License

Highlights from LINFACTOR: uses LU or CHOL to factorize a matrix, or previously computed factors to solve Ax=b

5.0
5.0 | 4 ratings Rate this file 12 Downloads (last 30 days) File Size: 6.96 KB File ID: #15138 Version: 1.2

LINFACTOR: uses LU or CHOL to factorize a matrix, or previously computed factors to solve Ax=b

Tim Davis (view profile)

29 May 2007 (Updated )

A simple M-file to solve Ax=b using LU or CHOL.

File Information
Description

LINFACTOR: factorize a matrix, or use the factors to solve Ax=b.
Uses LU or CHOL to factorize A, or uses a previously computed factorization to solve a linear system. This function automatically selects an LU or Cholesky factorization, depending on the matrix. A better method would be for you to select it yourself. Note that mldivide uses a faster method for detecting whether or not A is a candidate for sparse Cholesky factorization (see spsym in the CHOLMOD package, for example).

Example:

F = linfactor (A) ; % factorizes A into the object F

x = linfactor (F,b) ; % uses F to solve Ax=b

norm (A*x-b)

A second output is the time taken by the method, ignoring the overhead of determining which method to use. This makes for a fairer comparison between methods, since normally the user will know if the matrix is supposed to be symmetric positive definite or not, and whether or not the matrix is sparse. Also, the overhead here is much higher than mldivide or spsym.

This function has its limitations:

(1) determining whether or not the matrix is symmetric via nnz(A-A') is slow. mldivide (and spsym in CHOLMOD) do it much faster.

(2) MATLAB really needs a sparse linsolve. See cs_lsolve, cs_ltsolve, and cs_usolve in CSparse, for example. Using the sparse linsolve capability in KLU (part of SuiteSparse) can lead to a factor of 8 speedup in the forward/backsolve.

(3) this function really needs to be written as a mexFunction.

(4) the full power of mldivide is not brought to bear. For example, UMFPACK is not very fast for sparse tridiagonal matrices. It's about a factor of four slower than a specialized tridiagonal solver as used in mldivide.

(5) permuting a sparse vector or matrix is slower in MATLAB than it should be; a built-in linfactor would reduce this overhead.

(6) mldivide (x=A\b) when using UMFPACK uses relaxed partial pivoting and then iterative refinement. This leads to sparser LU factors, and typically accurate results. linfactor uses sparse LU without iterative refinement, and thus it can be slightly less accurate than x=A\b.

The primary purpose of this function is to answer The Perennially Asked Question (or The PAQ for short (*)): "Why not use x=inv(A)*b to solve Ax=b? How do I use LU or CHOL to solve Ax=b?" The full answer is below. The short answer to The PAQ (*) is "PAQ=LU ... ;-) ... never EVER use inv(A) to solve Ax=b."

The secondary purpose of this function is to provide a prototype for some of the functionality of a true MATLAB built-in linfactor function.

Finally, the third purpose of this function is that you might find it actually useful for production use, since its syntax is simpler than factorizing the matrix yourself and then using the factors to solve the system.

Oh, did I tell you never to use inv(A) to solve Ax=b?

Requires MATLAB 7.3 (R2006b) or later.

NOTE: If you have any difficulty with this file, or any of my files, please contact me directly rather than posting a "review". You cannot do any linear algebra on a cell array, to answer the one question below.

Acknowledgements

This file inspired Trilin and Don't Let That Inv Go Past Your Eyes; To Solve That System, Factorize!.

Required Products MATLAB
MATLAB release MATLAB 7.5 (R2007b)
MATLAB Search Path
/
/LINFACTOR
Other requirements License: this is free software with one condition: you must agree never to use inv(A) to solve Ax=b ... ;-)
12 Apr 2016 Gregory Vernon

10 Jun 2014 Amin

Amin (view profile)

11 Dec 2011 N. Yildirim

N. Yildirim (view profile)

14 Jun 2011 Katak Goreng

Katak Goreng (view profile)

i want to solve x=A\B where A is sparse, square,unsymmetric & have some zero diagonal whilst B is full vector. the size of A is >(100,000 x 100,000).what is the fastest solver to solve my problem? does matlab a\b is faster than KLU for this types of problem?

Comment only
04 Jun 2008 Ahmed Fasih

Note that this works when 'b' is a matrix as well! Works like a charm for Slepian-Bangs formulation for the Cramer-Rao bound!

07 Feb 2008 Shamil Akbar

>> F=linfactor(a)
??? Undefined function or method 'minus' for input arguments of type 'cell'.

Error in ==> linfactor at 118
if (nnz (A-A') == 0 & all (diag (A) > 0)) %#ok

HELP ME,,,,,, I WANT TO USE CHOLESKY TO FACTORIZE MY CORELATION MATRIX

Comment only
04 Jun 2007

Updated description; no change to code.

27 May 2009 1.1