Increasing accuracy of solving overdetermined linear system

Question

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?

Vector B: https://drive.google.com/open?id=18BFI9hd8d1MlwxjK9ByRSArhJ-h0rTWT
Matrix A (54x48): https://drive.google.com/open?id=1AaKEfz6fbPBbPutJSpSNlGou-ysLrAqr

Answer 1

matA = (* from your link *)
vecB = (* 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.