LU Factorization, Singularity Tolerance

I am computing an LU factorization with partial pivoting for dense matrices. I want to be able to detect when my matrix is singular to working precision. I am not sure what I should do for the tolerance on the diagonal values of U. What are your thoughts on what tolerance I should use?

For example

A = gallery('chebspec',3,0) % this will produce an ill-conditioned A. It has rank of 2.
[L,U,P] = lu(A)

A =

1.5 -2 0.5

0.5 0 -0.5

-0.5 2 -1.5

L =

1 0 0

-0.333333333333333 1 0

0.333333333333333 0.5 1

U =

1.5 -2 0.5

0 1.33333333333333 -1.33333333333333

0 0 -7.40148683083438e-17

P =

1 0 0

0 0 1

0 1 0

rank(A)

2

The U(3,3) is very small compared to the other diagonal values. It seems that something like a reasonable tolerance would be something like

d = abs(diag(U));
tolerance = 100*eps(max(d));
singularPivots = d < tolerance;

However, I think this idea is wrong compared to rank(). I need to either think this through or get help on what a good tolerance would be.

[Jason H. Nicholson]

Can I use the same tolerance that the MATLAB rank() for dense matrices uses? This seems reasonable. See code below:

A = gallery('chebspec',3,0); % this will produce an ill-conditioned A. It has rank of 2.
[L,U,P] = lu(A);
s = abs(diag(U));
tolerance = max(size(A)) * eps(max(s))
singularPivots = s < tolerance

tolerance =

6.6613e-16

singularPivots =

3×1 logical array

0

0

1

What do you think? Is this good enough? I am leaning towards this is good enough. My rational is eigen values of U are on the diaganol. The eigen value relative size to each other is similar but not the same as looking at the relative size of the singular values. However, the analog seems close enough that this may be good enough. Your thought are welcomed.

[Jason H. Nicholson]

This seems similar but slightly over my head:

A survey of condition number estimation for triangular matrices

[Jason H. Nicholson]

This seems to be a good resource: On the existence and computation of rank-revealing LU factorizations

I get the hint that the LU factorization I am getting from MATLAB may not reveal the rank. However, I just need approximate. Conservative is best meaning assuming it is rank deficient when maybe it isn't. Therefore, even if the tolerance I used in my previous comment was accurate, I am going to use 1000*tolerance to check for ill-conditioned and then take action with the bad rows/equations and small pivots/columns/variables accordingly.

[Jason H. Nicholson]

This is good too: LURP: Gaussian elimination with rook pivoting using efficient Fortran and MATLAB mex code

[Jason Nicholson]

Here is a good test matrix from Peters and Wilkinson:

n = 60;
A = eye(n) - triu(ones(n),1);
epsilon = -1/(2^(n-1)-1);
A(2:end,1) = epsilon;
rank(A) % 59
cond(A) % 3.6139e+18
[L,U,p] = lu(A,'vector');

s = abs(diag(U));
tolerance = max(size(A)) * eps(max(s)) % 1.3323e-14
singularPivots = s < tolerance % last pivot is singular
s(60) % 9.9920e-16

Peters, Gwen, and James Hardy Wilkinson. "The least squares problem and pseudo-inverses." The Computer Journal 13.3 (1970): 309-316.

[Jason Nicholson]

Another reasonable estimate of the degenerate pivot is

s = abs(diag(U)); % absolute value of pivots
tolerance = max(size(A))*norm(A,1)*eps(max(s));
singularPivots = s < tolerance

This one is based on code I found in null_lu.m in the Rank revealing lu decomposition submission:

tol         = max([n,m])*norm(A,1)*eps;

They just use eps. However, I think that isn't good enough. The eps value should be based on the largest pivot and thus eps(max(s));. This yields a larger tolerance for our degenerate test matrix.

n = 60;
A = eye(n) - triu(ones(n),1);
epsilon = -1/(2^(n-1)-1);
A(2:end,1) = epsilon;
rank(A) % 59
cond(A) % 3.6139e+18
[L,U,p] = lu(A,'vector');

s = abs(diag(U));
tolerance = max(size(A))*norm(A,1)* eps(max(s)) % 7.99360577730113e-13
singularPivots = s < tolerance % last pivot is singular
s(60) % 9.99200722162641e-16

[Jason Nicholson]

Four tests with our test matrix from Peters and Wilkinson:

• MATLAB lu() with partial pivoting. Idenifies the correct rank.
• MATLAB lu() with complete pivoting. Identifies the wrong rank.
• lurp() with p output. Identifies the wrong rank.
• lurp() with p and q output. Identifies the right rank.

test matrix:

n = 60;
A = eye(n) - triu(ones(n),1);
epsilon = -1/(2^(n-1)-1);
A(2:end,1) = epsilon;
rank(A) % 59
cond(A) % 3.6139e+18

MATLAB lu with partial pivoting

[L,U,p] = lu(A,'vector'); % partial pivoting
s = abs(diag(U));
tolerance = max(size(A))*norm(A,1)* eps(max(s)) % 7.99360577730113e-13
singularPivots = s < tolerance % last pivot is singular
s(60) % 9.99200722162641e-16

MATLAB lu with complete pivoting. This doesn't return the right singular pivots. Its not even close when comparing the tolerance and pivots. Note that the difference between the last singular value and and first 59 singular values is relatively large; this is not comparable to what the lu factorization pivots tell us for this case.

[L,U,p,q] = lu(sparse(A),'vector'); % force lu to use complete pivoting. Sadly, you have make A sparse to do this.
s = abs(full(diag(U)));
tolerance = max(size(A))*norm(A,1)* eps(max(s)) % 1.27897692436818e-11
singularPivots = s < tolerance % last pivot is singular
s(59:60) % [1.73472347597681e-18; 8.88178419700125e-16]
svd(A)

svd(A) =

37.2706744752901

12.4990715068938

7.58921991409077

5.51561758340263

4.38619501135786

3.68472771069481

3.21256290251896

2.87695954972578

2.62882661369057

2.4397752633566

2.29227608009941

2.1749454484734

2.08009011990592

2.00234024024841

1.93784631485687

1.88378610492719

1.83805042261933

1.79903655207152

1.76550874278132

1.73650179052771

1.71125303666314

1.68915354764701

1.66971250544675

1.65253086291163

1.63728160251772

1.62369477039393

1.61154600941432

1.60064768718541

1.59084196969222

1.58199536866755

1.57399441573992

1.56674220563485

1.56015561512363

1.55416305142467

1.54870261840903

1.54372061473283

1.53917029735375

1.53501085851238

1.53120657540464

1.52772610032638

1.52454186568362

1.5216295834053

1.51896782232114

1.5165376502368

1.51432232994807

1.51230706043397

1.5104787560675

1.50882585796973

1.5073381726731

1.50600673410593

1.5048236856013

1.50378217920187

1.50287629000192

1.50210094365879

1.50145185553412

1.50092548020233

1.50051897030098

1.50023014390266

1.50005745976793

1.03131216497561e-17

lurp() with p output. This is wrong. Not sure what to say here. There are not 51 singular pivots.

[L,U,p] = lurp(A,'vector');

s = abs(diag(U));
tolerance = max(size(A))*norm(A,1)* eps(max(s)) % 7.99360577730113e-13
singularPivots = s < tolerance % pivots 10:60 are considered singular
s

s =

1

1

1

1

1

1

1

1

1

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

0

lurp with p and q outputs. This gets it right.

[L,U,p,q] = lurp(A,'vector');

s = abs(diag(U));
tolerance = max(size(A))*norm(A,1)* eps(max(s)) % 1.59872115546023e-12
singularPivots = s < tolerance % last pivot is singular
s(60) % 9.95936892841527e-30