For the case of PA=LU factorization, I found some documents which tell that it may delete the probability of having 0's on the diagonal of Matrix A. But I am not sure if I got it right. If so, what is the problem of having 0's on the diagonal of Matrix A?

Other questions which are related, to which I couldn't find any easy explanation:

  • What are the advantages / disadvantages of LU factorization? And when to use it?
  • What are the advantages / disadvantages of PA=LU factorization? And when to use it?
  • What are the advantages / disadvantages of Cholesky factorization? And when to use it?
  • What are the advantages / disadvantages of QR factorization? And when to use it?

Thank you for any help!

  • 1
    $\begingroup$ A good textbook on numerical linear algebra would contain the answers to all your questions. There is a reference-request tag which you can use to ask for recommendations. $\endgroup$ – Carl Christian Jun 30 '18 at 16:29
  • What are the advantages / disadvantages of LU factorization? And when to use it?
  • In practice for when you use the solver on your computer it is going to go through a process and determine whether it can first use Cholesky then LU to solve the system of equations. In general the normal LU decomp isn't used because of the stability.
  • What are the advantages / disadvantages of PA=LU factorization? And when to use it?
  • Some advantages are that it is more numerically stable and it also tells you the rank of the matrix.
  • What are the advantages / disadvantages of Cholesky factorization? And when to use it?
  • You'd use it if your matrix was positive definite. You can't otherwise. It gives a speed up.
  • What are the advantages / disadvantages of QR factorization? And when to use it?
  • You would utilize the QR decomposition for numerical stability. It is the back bone of most techniques for orthogonalization. There are many QR decomps but you're really referring to Householder. LU is not backwards stable.
  • When would use SVD ?
  • The SVD would be used to control the condition number in a decomposition. It is only marginally more time consuming than the QR.

These are some examples here...ok. For the LU decomposition there exists matrices which are pathologically ill-conditioned. There is a book on them by Nicholas Higham. He published a matrix testing thing to demonstrate some of this, however some of it still exists in Matlab. I can't find these at the moment. It's particularly hard to demonstrate this as the solver are so good. Matlab has a chart which follows like. enter image description here

Simply, for the questions.

The PLU is more numerically stable than the general LU decomposition. It actually isn't used in practice. It is what they are referred to. Following here.

function [L,U,P] = my_lu_piv(A)
n = size(A,1);
I = eye(n);
O = zeros(n);
L = I;
U = O;
P = I;
    function change_rows(k,p)
        x = P(k,:); P(k,:) = P(p,:); P(p,:) = x;
        x = A(k,:); A(k,:) = A(p,:); A(p,:) = x;
        x = v(k); v(k) = v(p); v(p) = x;

    function change_L(k,p)
        x = L(k,1:k-1); L(k,1:k-1) = L(p,1:k-1);
        L(p,1:k-1) = x;
    for k = 1:n
        if k == 1, v(k:n) = A(k:n,k);
            z = L(1:k-1,1:k -1)\ A(1:k-1,k);
            U(1:k-1,k) = z;
            v(k:n) = A(k:n,k)-L(k:n,1:k-1)*z;
        if k<n
            x = v(k:n); p = (k-1)+find(abs(x) == max(abs(x))); % find index p
            L(k+1:n,k) = v(k+1:n)/v(k);
            if k > 1, change_L(k,p); end
        U(k,k) = v(k);

function [L,U] = my_LU(A)
n = size(A,1);
I = eye(n);
L = I;
U = A;
    for k=1:n-1

        L(k+1:n,k) = U(k+1:n,k)/U(k,k);
        for j=1:k+1:n
            U(j,k:n) = U(j,k:n)-L(j,k)*U(k,k:n);


    A = [2,1,1,0;4,3,3,1;8,7,9,5;6,7,9,8];
    [L,U,P] = my_lu_piv(A);
    [L1,U1] = my_LU(A);
    c = cond(A)
    err1 = norm(A-P*L*U);
    err2 = norm(A-L1*U1);
   c =


   err1 =


   err2 =


The advantage of the knowing the rank of the matrix generally people want to know the rank of the matrix. There is a command for the rank. However, I believe that newer techniques are based on low-rank approximations. E.g. like the following.

 function A = genrankk(m,k,n)
         % Generates rank k matrix
         A = rand(m,k)*rand(k,n);


     A = genrankk(m,k,n);

    [U,S,V] = svd(A);
    [m,n] = size(S,'econ');
    S =

    5.0708         0         0         0         0         0         0
         0    0.4292         0         0         0         0         0
         0         0    0.1937         0         0         0         0
         0         0         0    0.0000         0         0         0
         0         0         0         0    0.0000         0         0



Now these are machine precision because the rank is actually 3. To generate a low rank rank approximation we can do the following..I technically did but..

function B = makeapprox(U,S,V,k)
% rank k approx

B = U(:,1:k)*S(1:k,1:k)*V(:,1:k)';


   function [l1,l2] =  errf(A)
    % Vector of error between approximation

    [m,n] = size(A);
    [U,S,V] = svd(A);
    l = min(m,n);
    l1 = linspace(1,l,l);
    l2 = l1;
    for i =1:l
        B = makeapprox(U,S,V,i);
        l1(i) = norm(A-B);


 [l1,l2]  = errf(A)

Technically the QR decomposition should be used when the condition number of the matrix is high. E.g. you can generate a matrix with a bad condition number

function A = randommatrix(m,n,kappa)
% Returns an m x n matrix with condition number kappa,
% exponentially graded singular values
k = min(m,n);
[Q1,R1] = qr(randn(m,m));
[Q2,R2] = qr(randn(n,n));
% Originally thought we would evaluate for any m x n matrix, so if not square took minimum dimension
U = Q1(1:m,1:k);
V = Q2(1:k,1:n);
% Take k-1 root of condition number to get number for grade
l = kappa^(1/j);
% Construct the singular value matrix, l^0 = 1, then l^(k-1)= 1/kappa 
S = diag(l.^(0:-1:-j));

A = U*S*V;


and the QR decomp should perform better at solving the system

  • 1
    $\begingroup$ Technically PA=LU is actually backward stable, just the constant factor is highly non-uniform in the dimension. $\endgroup$ – Ian Jun 30 '18 at 13:13
  • $\begingroup$ @Metahominid thanks for you answer! 1) You say PA=LU is more numerically stable. More than what? Than LU? 2) What would be a "concrete" example of unstability of LU vs Choles.? 3) As of your answer to PA=LU, what is the advantage of knowing the rank of the Matrix? 4) What is the advantage of having a positive definite matrix? Could you give a "concrete" example, of a matrix, where you would rather use QR decomp? I get a bit the general picture, but can't tell in a concrete example when to use which decomposition! It would be REALLY helpful, if you could provide concrete examples when use what? $\endgroup$ – ZelelB Jul 2 '18 at 18:44
  • 1
    $\begingroup$ @Zeleb 1. More numerically stable than standard LU decomp. One second I'll update it with some examples. I don't you'd have a trade off vs. LU and cholesky in stability. The trade off is in speed. Cholesky is about twice as fast as the standard LU decomp. They are the same basic decomposition. 3. Typically low rank approximations. 4. Yes, one moment. $\endgroup$ – воитель Jul 2 '18 at 18:49
  • $\begingroup$ @Metahominid Man! Thank you so much for this great answer and the top explanation! Really! $\endgroup$ – ZelelB Jul 5 '18 at 9:55
  • $\begingroup$ @Metahominid The graph is so great! Just three small things: 1) After "Is A upper Hessenberg?" --> No, wouldn't you check if "A is well-conditioned?" if "Yes", LU-decompostion, if "No", QR (with Householder)? 2) for "Are all diagonal elements positive or all negative?" before doing Cholesky: I thought we just check if "All diagonal elements are "POSITIVE". If all are negative, could we still do Cholesky? 3) What do you mean by "Does Cholesky succeed?" How to check that? $\endgroup$ – ZelelB Jul 5 '18 at 10:03

Partial answer:

If you have a system $A x = b$, you can decompose the matrix $A$ to find $x$ that minimizes $\|Ax-b\|_2$. Analytically, the solution is $x = A^\dagger b$ where $A^\dagger$ is the pseudo-inverse of $A$. However, determining the pseudoinverse of $A$ analytically (perhaps as $A^\dagger = (A^T A)^{-1} A^T$ is numerically unstable.

So you can set $A=LU$ or $A=QR$ or $A=U\Sigma V^T$ where the last decomposition is the singular value decomposition. These decompositions permit numerically stable estimations of the psuedo-inverse, and can be used to solve the optimization problem. $LU$ is the fastest decomposition, but the least numerically stable. The SVD is the slowest decomposition, but permits the most numerically stable algorithm.

  • $\begingroup$ thank you for your answer! Could you please provide an example, where you would use QR rather than LU and another where you would use LU instead of QR? That would really help me, because in theory it is quite clear! But in exercices, I couldn't really tell, when to use what, and in which case would be which decomposition better and why. Thanks again! $\endgroup$ – ZelelB Jul 2 '18 at 18:46
  • 1
    $\begingroup$ @ZelelB In my mind, it depends on your application. Are you working in medicine or space exploration where a failure of factorization will lead to death? Then you probably want to use the SVD. Are you working on a new phone app game and you need something quick that works on minimal hardware? Then LU seems like a good way to go! In "Numerical Linear Algebra", Trefethen recommends QR for "daily use". $\endgroup$ – NicNic8 Jul 3 '18 at 1:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.