Ggtkuk

I have some composition data for a geological sample and its constituent minerals. I want to estimate the proportions of the minerals that make up the bulk sample. This is essentially a mass balance problem. Each variable is pseudo-independent (but constrained by constant sum).

I'm starting with the matrix form $a.x=b$,
where, $a$ is the array $elemental~composition~vector times mineral ~species$, $x$ is the column vector of mineral proportions, and $b$ is the bulk composition. The key constraint with this kind of problem is that

$0 leq x_j leq 1~~~~~(j=1,...m)$

and

$sum_j=1^m x_j=1$

here is some example data

a=63.1545, 64.3049, 100., 37.1417, 32.4026, 30.0382, 30.7033, 0., 36.5444, 0., 0., 0., 0.0034,
 0.0016, 0.01, 0., 0.6641, 2.35946, 26.91, 0., 0., 0.297125, 0., 0., 0., 51.3026, 
 17.8683, 21.096, 0., 11.2387, 14.3019, 6.46115, 0.34805, 0., 14.368, 0., 0., 0., 0., 
 0.0223, 0.1191, 0., 30.6859, 32.0779, 1.5471, 1.5014, 0., 11.5299, 0., 0., 0., 43.3621, 
 0., 0., 0., 0.94408, 0.6131, 0.0986, 0., 0., 0., 0., 0., 0., 2.1039, 
 0., 0., 0., 1.14598, 1.82486, 0.03695, 0.00935, 0., 0.2312, 0., 0., 0., 0.0166, 
 0.625667, 9.77073, 0., 1.35584, 0.0558, 0.05635, 0., 0., 0.079275, 0., 0., 0., 0.018, 
 0., 2.80603, 0., 10.7388, 0.0245429, 27.7454, 2.19443, 55.08, 10.5478, 56.03, 78.13, 10.44, 0.0262, 
 15.7365, 0.0874667, 0., 1.89652, 8.54851, 0.00475, 0., 0., 0.09665, 0., 0., 0., 0.0115, 
 0., 0., 0., 0.16636, 0., 0., 0.0293, 42.4, 0.0114, 0., 0., 0., 0.0617

b=65.67, 0.52, 14.77, 5.418, 0.13, 0.3, 3.22, 2.05, 6.07, 0.13

If I use the LeastSquares function I get

lsq = LeastSquares[a, b]

0.331681, 0.283916, 0.166439, 0.172393, 0.0586919,
0.0795275,
0.00802007, 0.00244529, -0.030593, -0.0157832, -0.0220087,
-0.00294087, -0.0363984

This result looks OK, except negative proportions are not allowed and the sum does not equal 1.

Total@lsq

0.99539

How do I implement linear least squares fitting for multivariate data under the constraint that values lie between 0 to 1 and sum to 1 in Mathematica?

Also how can I derive the residuals for the fit?

One approach is to reframe it as a minimization problem:

xVec = Array[x, 13];
NMinimize[Total[(a.xVec - b)^2], Total[xVec] == 1, Thread[xVec >= 0], xVec]

You can also get an improved residual by relaxing the equality constraint:

NMinimize[Total[(a.xVec - b)^2] + (Total[xVec] - 1)^2, Thread[xVec >= 0], xVec]

Thanks to Roman for some simplifications.

This can be done as a linear programming problem although the optimization is not least-squares in that case. The idea is to set up new variables, which are absolute values of discrepancies between a.x. and b. This can be done as below.

xvec = Array[x, Length[a[[1]]]];
dvec = Array[d, Length[a]];
lpolys = a.xvec - b;
ineqs = Flatten[Thread[xvec >= 0], Thread[dvec - lpolys >= 0], 
 Thread[dvec + lpolys >= 0], Total[xvec] == 1];

Now use NMinimize.

min, vals = 
 NMinimize[Total[dvec], ineqs, Join[xvec, dvec], PrecisionGoal -> 10]
Total[xvec /. vals]

(* Out[90]= 0.627674050324, x[1] -> 0.339337833732, 
 x[2] -> 0.292918812222, x[3] -> 0.195908896153, 
 x[4] -> 0.10491683141, x[5] -> 0.0591148361052, x[6] -> 0., 
 x[7] -> 0., x[8] -> 0.00181140763364, x[9] -> 0., x[10] -> 0., 
 x[11] -> 0., x[12] -> 0., x[13] -> 0.00599138274489, d[1] -> 0., 
 d[2] -> 0., d[3] -> 0.502611255236, d[4] -> 0., 
 d[5] -> 0.0178984583702, d[6] -> 0.0717916527728, d[7] -> 0., 
 d[8] -> 0., d[9] -> 0., d[10] -> 0.0353726839451

Out[91]= 1. *)