# Choosing appropriate WorkingPrecision when solving numerical system of equations [closed]

Consider a linear system of equations, which is conveniently written as A.x=y. The matrix A has dimensions 500x500, while x and y are of course 500 dimensional vectors. Matrix A and vector y contain numerical (float) numbers that in principle can be generated with arbitrary precision. Analytically, it is expected that x will contain mostly zeros and a couple of integers of magnitude <10. It turns out that a generic matrix A, in the particular problem at hand, has condition number such that:

Log10[LinearAlgebraMatrixConditionNumber[A]]


27.053931

which means that while x is being determined, around 27 digits of precision are expected to be lost purely due to the properties of matrix A (even before the loss due to rounding errors etc). Considering the above, what starting precision for entries of matrix A and vector y should I shoot for in order to unambiguously determine vector x while spending as little calculation time as possible? Is there any estimate one can make apart from doing several experimental trial and error calculations?

EDIT

I have been asked where the matrix A comes from. It's origin is in Physics. Here is how it is constructed. First we generate n=6 light-like Minkowski-space momenta which add up to zero:

n = 6; precis = 200; max = 1;
Do[
Subscript[p, i, 1] = RandomReal[{-max, max}, WorkingPrecision -> precis];
Subscript[p, i, 2] = RandomReal[{-max, max}, WorkingPrecision -> precis];
Subscript[p, i, 3] = RandomReal[{-max, max}, WorkingPrecision -> precis];
Subscript[p, i, 0] = (-1)^RandomInteger[{0, 1}] Sqrt[Subscript[p, i, 1]^2 + Subscript[p, i, 2]^2 + Subscript[p, i, 3]^2];
, {i, 1, n - 2}];
Subscript[p, n - 1, 2] = -Sum[Subscript[p, i, 2], {i, 1, n - 2}];
Subscript[p, n - 1, 3] = -Sum[Subscript[p, i, 3], {i, 1, n - 2}];
Subscript[p, n, 2] = 0; Subscript[p, n, 3] = 0;
AA = Sum[Subscript[p, i, 0], {i, 1, n - 2}];
BB = Sum[Subscript[p, i, 1], {i, 1, n - 2}];
DD = Subscript[p, n - 1, 2]^2 + Subscript[p, n - 1, 3]^2;
Subscript[p, n, 0] = (-AA^2 + BB^2 + DD)/(2 (AA + BB));
Subscript[p, n, 1] = -((-AA^2 + BB^2 + DD)/(2 (AA + BB)));
Subscript[p, n - 1, 0] = -(((AA + BB)^2 + DD)/(2 (AA + BB)));
Subscript[p, n - 1, 1] = (-(AA + BB)^2 + DD)/(2 (AA + BB));
Table[Sum[Subscript[p, i, j], {i, 1, n}], {j, 0, 3}]
Table[Subscript[p, i, 0]^2 - (Subscript[p, i, 1]^2 + Subscript[p, i, 2]^2 + Subscript[p, i, 3]^2), {i, 1, n}]
pq[x_, y_] := Subscript[p, x, 0] Subscript[p, y, 0] - (Subscript[p, x, 1] Subscript[p, y, 1] + Subscript[p, x, 2] Subscript[p, y, 2] + Subscript[p, x, 3] Subscript[p, y, 3]);


As the output from the two tables demonstrates, the light-like and add-up-to-zero constraints are properly satisfied. Function pq[x,y] gives the inner product over Minkowski space. Next, define a function that when called will change the momenta slightly but maintain the conditions:

ChangeMomentaSlightly[] := Module[{},
max = 1/10000;
Do[
Subscript[p, i, 1] = RandomReal[{Subscript[p, i, 1] - max, Subscript[p, i, 1] + max}, WorkingPrecision -> precis];
Subscript[p, i, 2] = RandomReal[{Subscript[p, i, 2] - max, Subscript[p, i, 2] + max}, WorkingPrecision -> precis];
Subscript[p, i, 3] = RandomReal[{Subscript[p, i, 3] - max, Subscript[p, i, 3] + max}, WorkingPrecision -> precis];
Subscript[p, i, 0] = Sign[Subscript[p, i, 0]] Sqrt[Subscript[p, i, 1]^2 + Subscript[p, i, 2]^2 + Subscript[p, i, 3]^2];
, {i, 1, n - 2}];
Subscript[p, n - 1, 2] = -Sum[Subscript[p, i, 2], {i, 1, n - 2}];
Subscript[p, n - 1, 3] = -Sum[Subscript[p, i, 3], {i, 1, n - 2}];
Subscript[p, n, 2] = 0; Subscript[p, n, 3] = 0;
AA = Sum[Subscript[p, i, 0], {i, 1, n - 2}];
BB = Sum[Subscript[p, i, 1], {i, 1, n - 2}];
DD = Subscript[p, n - 1, 2]^2 + Subscript[p, n - 1, 3]^2;
Subscript[p, n, 0] = (-AA^2 + BB^2 + DD)/(2 (AA + BB));
Subscript[p, n, 1] = -((-AA^2 + BB^2 + DD)/(2 (AA + BB)));
Subscript[p, n - 1, 0] = -(((AA + BB)^2 + DD)/(2 (AA + BB)));
Subscript[p, n - 1, 1] = (-(AA + BB)^2 + DD)/(2 (AA + BB));]


Now, using s[Subscript[k, x] + Subscript[k, y]] as an analytic representation for pq[x,y] (and equivalently pq[y,x]) we can construct one particular set of linearly independent inner products among the n vectors (this gives 9 different elements):

mySet = DeleteCases[ DeleteCases[Table[s[Subscript[k, x] + Subscript[k, y]], {x, 2, n}, {y, 2,n}] /. s[2 Subscript[k, _]] -> 0 // Flatten //DeleteDuplicates, 0], s[Subscript[k, 2] + Subscript[k, 3]]];


Then create a table of all possible unique multiplications of four elements in mySet with each other (this gives 495 different elements):

myTab = Table[ mySet[[i]] mySet[[j]] mySet[[q]] mySet[[w]], {i, 1, 9}, {j, 1, 9}, {q, 1, 9}, {w, 1, 9}] // Flatten // DeleteDuplicates;


For the sake of discussion, let's create a sample vector y which we can experiment with:

rndm[x_, y_] := RandomInteger[{x, y}];
sampleY = Sum[rndm[0, 1] rndm[1, 10] myTab[[rndm[1, 495]]], {i, 1, rndm[1, 495]}];


Finally, numerical matrix A is created as:

A = DiagonalMatrix[Range[495]];
y = Table[0, {i, 1, 495}];
Do[ChangeMomentaSlightly[];
y[[i]] = sampleY /. s[Subscript[k, x_] + Subscript[k, y_]] :> pq[x, y];
Do[A[[i, j]] = myTab[[j]] /. s[Subscript[k, x_] + Subscript[k, y_]] :> pq[x, y];, {j, 1, 495}];
, {i, 1, 495}]


The actual vector y is not known analytically and there is a certain lengthy calculation that can be used to obtain it numerically. The whole point of the above procedure is to fix the coefficients in front of a complete set of functions in order to determine what y is analytically. The sample data above seems to work out fine at a working precision of 200:

fun = LinearSolve[A];
fun[y] // Chop // Rationalize


However, it would be nice to be able to estimate the sweet spot for working precision in order to automate the process and spend as little computation time as possible on each A and y`.

## closed as off-topic by MarcoB, m_goldberg, user9660, Öskå, JensMar 31 '16 at 15:43

• The question does not concern the technical computing software Mathematica by Wolfram Research. Please see the help center to find out about the topics that can be asked here.
If this question can be reworded to fit the rules in the help center, please edit the question.

• Apart from sparsity, is there any other structure in your matrix? Would you mind telling where it came from? – J. M. will be back soon Aug 5 '15 at 7:15
• I'm voting to close this question as off-topic because it is too localized; i.e, it applies only to the local situation and needs of its poster and answers will not benefit others. – m_goldberg Mar 29 '16 at 14:11