# Determining multivariate least squares with constraint

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?

• My usual soapbox speech: Why not consult a statistician first and then figure out how to implement an appropriate process in Mathematica? Consider stats.stackexchange.com/questions/267014/…. – JimB Mar 20 at 14:53
• @JimB precisely because someone will tell me that I can't use continuum methods with compositional data. The constant sum problem is considered a vague oddity in geological circles and largely ignored (for better or worse). – geordie Mar 20 at 15:26
• I'm sure there must be more to it than "They won't let me do what I want!" Compositional data analysis has a pretty standard set of statistical approaches. Chemists use it all the time. These approaches are many times called "mixture designs": itl.nist.gov/div898/handbook/pri/section5/pri54.htm. Maybe you've been talking to the wrong folks. – JimB Mar 20 at 15:49
• Like I said, geologists don't. Perhaps it has to do with the general 'noisiness' of geoscience datasets. Most other sciences are able to rely on fairly precise numbers. – geordie Mar 20 at 16:22
• "Most other sciences are able to rely on fairly precise numbers." Ha, ha. You're funny! – JimB Mar 20 at 16:30

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.

• Thread[xVec > 0] would be simpler on the third line and does the same thing. Also, you can constrain the sum by equality: NMinimize[{(a.xVec-b).(a.xVec-b), Total[xVec]==1, Thread[xVec>0]}, xVec]. I think the square in Total[(a.xVec - b)]^2 should be taken inside of the Total, not outside. – Roman Mar 20 at 14:15

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;
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. *)