Numerical Error with Large Matrices

Question

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?

Answer 1

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[[1]]

(* Out[275]= {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[278]:= prod3 = inv2.matrix; Max[Abs[prod3[[1]]]]

Out[279]= 55.8571828306

That's a far cry from the value 1. And indeed the condition number is not good.

LUDecomposition[mat2][[3]]

(* Out[280]= 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[283]= {1.27696286402*10^21} *)

SingularValueList[mat2, -1]

(* Out[285]= {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.

Answer 2

this worked for me

    size = 100; 
    min = 100000000; 
max = 1000000000; 
ltm = Table[If[j <= i, RandomReal[
      {min, max}], 0], {i, size}, 
    {j, size}]; 
matrix = Rationalize[
    ltm . Transpose[ltm], 
    10^(-30)]; 
MatrixForm[inv = N[Inverse[matrix]]]
Chop[inv, 10^(-10)]

Hoping this approach helps you.