# Increasing accuracy of solving overdetermined linear system

I am given $48 \times 48$ matrix $A$ and a vector $b$ and I would like to solve system $Ax = b$. I know that $A$ is underdetermined, i.e. there exist many solutions for $x$. Due to some considerations, I also know that certainly, six additional equations must hold (by adding which system of equations becomes unique).

Roughly, I am calculating probabilities where $x$ contains probabilities for each of 48 events. These six additional equations are because some of the groups of events are independent and equally likely, so I know that each of this group has probability $1/6$.

To implement this in Mathematica I added 6 more rows to $A$, where I put $1$ if the corresponding element in $x$ is a member of some group and $0$ for all that are not members of some group. Then, for $b$ I added 6 numbers, where each of them is $1/6$.

Therefore, I have an overdetermined equation, which I know must have a unique solution. To solve it, I used LeastSquares method. To check whether the solution is accurate enough, I computed the total probability of all events (which by 6 constraint equations should be equal to 1). I found out that computation leads to the value $0.53$, which is unsatisfactory for me. Then, I tried minimizing 1-norm and that leads to the value $0.84$.

My question:

1) What is the most efficient way to solve this problem? I understand that numerical calculations are never perfect, but is there a way to tell Mathematica to consider constraint equations as precise as possible and let numerical errors go into original 48 equations? Having total probability close to $1$ is crucial to me.

2) If I know that my system of overdetermined equations is unique, is there any way of solving it precisely using Mathematica? I.e., what is the way to reduce numerical errors as much as possible?

• Hard to answer without concrete matrices. At 1) Sounds as if you were looking for LeastSquares. At 2) LinearSolve with exact input will try to compute exact solutions. With floating point input and such tiny matrices, it is quite unlikely that the numerical errors will be be measurable. – Henrik Schumacher Jul 28 '18 at 16:06
• @HenrikSchumacher Thanks for the comments, I am not sure I understand the 2). LinearSolve (as I understand) solves only for quadratic matrices. When I add my constraints then the matrix is $54 \times 48$, and LinearSolve does not give answer. – Daniels Krimans Jul 28 '18 at 16:22
• Ah, I read "unique" for solvable. Anyways, LeastSquares should do what you want. – Henrik Schumacher Jul 28 '18 at 16:31
• If the right hand side lies within the range of A, LeastSquares returns the solution with the least $2$-norm. So, if the system is solvable and if there is only one solution then both LinearSolve and LeastSquares will return the unique solution. – Henrik Schumacher Jul 28 '18 at 18:33
• (I was starting to write substantially the same as @HenrikSchumacher, so I'll just add...) This does not strike me as something one should be embarassed about. A least-squares solution can be useful in its own right even when a system is inconsistent. What you might want to check is whether you got the model slightly wrong, or have a coding error in generating the matrix, etc. – Daniel Lichtblau Jul 28 '18 at 18:37

matA = (* from your link *)
xVec = ToExpression[Table["x" <> ToString[i], {i, 48}]];


Choose six variables such that the corresponding columns in the matrix of constraint-rows has full rank:

elimVars = {3, 4, 5, 6, 7, 8};
If[MatrixRank[matA[[-6 ;;, elimVars]]] < 6, "CHOOSE OTHER ELIMVARS"]
(* Null *)


Eliminate the variables and substitute the solution in matA and vecB

elim = Solve[matA[[-6 ;;]].xVec == vecB[[-6 ;;]], xVec[[elimVars]]];

newMat = Transpose[With[{newExpr = matA.(xVec /. First[elim])},
newb = vecB - newExpr /. Thread[xVec -> 0];
Coefficient[newExpr, #] & /@ Delete[xVec, Transpose[{elimVars}]]]];

res = xVec /. First[elim] /. Thread[Delete[xVec, Transpose[{elimVars}]] ->
LeastSquares[newMat, newb]];

matA[[-6 ;;]].res


{0.16666667, 0.16666667, 0.16666667, 0.16666667, 0.16666667, 0.16666667}

res contains the probabilities, and it is the same however elimVars is chosen as long as the submatrix has full rank.

Alternative:
The modification below constrains only the sum of the event's probabilities. The result is very similar:

xVec = ToExpression[Table["x" <> ToString[i], {i, 48}]];
elimVar = 1;
elim = Solve[Total[matA[[-6 ;;]].xVec] == 1, xVec[[elimVar]]];
newMat = Transpose[With[{newExpr = matA.(xVec /. First[elim])},
newb = vecB - newExpr /. Thread[xVec -> 0];
Coefficient[newExpr, #] & /@ Delete[xVec, elimVar]]];

res = xVec /. First[elim] /. Thread[Delete[xVec, elimVar] -> LeastSquares[newMat, newb]];

matA[[-6 ;;]].res
Total[matA[[-6 ;;]].res]


{0.16671407, 0.16755255, 0.16477448, 0.16694582, 0.1663015, 0.16771158}
1.

• Thanks for the answer. I am still going through it. Just a fast question - is it possible to solve something similar if groups are given by some equation? For example, if I have two groups and x = {x1, x2}, then 0.4*x1+0.6*x2 = 1/2 and 0.5*x1+0.8*x2 = 1/2. In a sense, groups might overlap and group sums are conditions with coefficients. – Daniels Krimans Jul 28 '18 at 16:45