Numerical Error with Large Matrices

I am writing a Finite Element Analysis program in Mathematica. The code involves handling a large matrix with large entries. I get an error when I try to use Mathematica's "LinearSolve" to solve a problem that requires numbers larger than approximately 1e7 in the matrix; the error states that I have a "badly conditioned matrix".

As a test, I generated an arbitrarily sized lower triangular matrix with random entries between a prescribed min and max value. Then I form an invertible matrix by multiplying the lower triangular matrix by its transpose (it can be shown that a lower triangular matrix multiplied by its transpose is invertible). Then, I multiply the newly formed matrix by its inverse. The result should be the identity matrix. Just as I suspected, when the matrix is either large or contains large entries (i.e. 50x50 with numbers between 1e7 and 1e8), I get the "badly conditioned matrix", and "may contain numerical errors" warnings.

Here is the code:

ClearAll["Global*"];
size = 50;
min = 10*^7;
max = 10*^8;
ltm = Array[0 &, {size, size}];
matrix = Array[0 &, {size, size}];
For[i = 1, i <= size, i++,
For[j = 1, j <= i, j++,
ltm[[i, j]] = RandomReal[{min, max}];
]
]
matrix = ltm.Transpose[ltm];
inv = Inverse[matrix];
Print[matrix.inv // MatrixForm];

I can fix the error by forcing the random entries to have a set precision of 30:

ClearAll["Global*"];
size = 50;
min = 10*^7;
max = 10*^8;
ltm = Array[0 &, {size, size}];
matrix = Array[0 &, {size, size}];
For[i = 1, i <= size, i++,
For[j = 1, j <= i, j++,
ltm[[i, j]] = SetPrecision[RandomReal[{min, max}],30];
]
]
matrix = ltm.Transpose[ltm];
inv = Inverse[matrix];
Print[matrix.inv // MatrixForm];

The problem is that there is no easy place in my actual program to use the SetPrecision[] command. I tried to use

$MinPrecision=30; at the top of the sample code, but it does not work. Any thoughts? 1 Answer Can do as follows. First generate the matrix of machine numbers. size = 100; min = 10*^7; max = 10*^8; ltm = Table[ If[j <= i, RandomReal[{min, max}], 0], {i, size}, {j, size}]; matrix = ltm.Transpose[ltm]; Now set the precision, invert, check the product. use the higher prec matrices throughout. mat2 = SetPrecision[matrix, 30]; inv2 = Inverse[mat2]; prod2 = inv2.mat2; prod2[] (* Out= {1.000000000000, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13, 0.*10^-13} *) To give an idea of the error amplification (which the condition number gauges), instead multiply this carefully computed inverse by the original matrix. In:= prod3 = inv2.matrix; Max[Abs[prod3[]]] Out= 55.8571828306 That's a far cry from the value 1. And indeed the condition number is not good. LUDecomposition[mat2][] (* Out= 2.89055384099761752993388785181*10^21 *) We'll also check the condition number in the 2-norm. It is the quotient of the largest and smallest singular values. SingularValueList[matrix, 1] (* Out= {1.27696286402*10^21} *) SingularValueList[mat2, -1] (* Out= {0.825206647809554673903655166033} *) Note that in order to compute this smallest one I needed to use the higher precision matrix. The bottom line (located conveniently at the bottom) is that you will need to be quite careful in how you use these matrices because the error sensitivity will be considerable. About the last thing I would recommend is actually inverting them. No,actually that's next-to-last; taking the determinant and hoping to get a valid result would be last. • I wrapped the final coefficient matrix with SetPrecision[] before using it in the LinearSolve[] function and that solved my problem. My question is why does Mathematica not automatically use that precision, or why does the error occur when I define the$MinPrecision to be 30 (even though I'm pretty sure it is defined higher by default)? – Matthew Apr 3 '13 at 13:18
• (1) Mathematica will not change the precision of your input unless directed to do so e.g. with SetPrecision. In particular, operations involving machine doubles will almost always be computed using same. – Daniel Lichtblau Apr 3 '13 at 15:22
• (2) The problem is not necessarily "solved" in the way you claim. What you did is create (via Setprecision a particular matrix in a continuous range of matrices, and do a computation with that matrix. It is not necessarily going to correspond well with analogous computations using different matrices in that same initial range. If your input is only known to a few decimal places then you might or might not be computing a viable result. – Daniel Lichtblau Apr 3 '13 at 15:26