matrix math: finding random solution

Share icon

Wed, 06/19/2013 - 05:30 am

I have a NxM matrix AA and a wave bb with M rows. Now I'd like to find random elements in a wave cc with N rows such that, if I execute:

function foo (AA, cc)
    wave AA, cc
   
    Duplicate/FREE/O AA, frac
    frac = AA[p][q] *cc[p]
    MatrixOP/O bb_calc = (SumCols(frac))^t 
end



bb_calc and bb are identical. A constraint is that all elements in cc must be > 0.

E.g., if AA is:
aa[0][0]= {1,2,2,0,0,1,0,0,0,1}
aa[0][1]= {0,0,0,1,2,1,1,0,0,0}
aa[0][2]= {0,0,1,0,0,0,1,1,2,1}

and bb:
bb[0]= {2,1,1}

a possible solution would be:

cc[0]= {0.5,0.15,0.5,0.5,0.15,0.1,0.1,0.05,0.125,0.1}


I'd be grateful about ideas and suggestions!
Hello Charlie,

Your function foo() looks like simple matrix multiplication of the transpose of AA so you can replace it with:

MatrixOP/O bb_calc = (AA^t) x cc


Now that you cast it into this form, it appears that you are looking for a solution of a system of the form:

MB=MA x CC

with MA (3x10 matrix) and CC (10,1).

If you were to use MatrixLLS on MA and MB you would get a minimum norm solution -- not one that restricts all elements of CC to be positive.

Clearly you could throw this into some optimize or root finding code but I'd be curious if there are other properties of the matrix AA that one could use. In particular, if AA is not random, I'd start by looking at its SVD.

AG

June 19, 2013 at 11:02 am - Permalink

Igor wrote:
Hello Charlie,

Your function foo() looks like simple matrix multiplication of the transpose of AA so you can replace it with:

MatrixOP/O bb_calc = (AA^t) x cc



Hi AG,
yes, of course! I ignored that there is "x" and "*" in MatrixOP! D'oh! I thought about MatrixLLS, but that doesn't solve the non-negative, non-zero issue.


Quote:

Clearly you could throw this into some optimize or root finding code but I'd be curious if there are other properties of the matrix AA that one could use. In particular, if AA is not random, I'd start by looking at its SVD.


This is the background to the question: I want to calculate chemical equilibrium using a Gibbs free energy minimisation technique, which however requires an initial guess for the abundance of the different species in the system. For simple systems this is easily done by hand, but of course an auto-guess is desirable. This guess must satisfy the mass balance constraint. AA is in this case a composition matrix with species in rows (e.g. CO, CO2, H2O, CH4, O2 ...) and number of atoms per species in columns (e.g. C, O, H), and bb is the total number of C, O, H in the system. cc represents the moles of each species.

Is there any way to tweak MatrixLLS, such that it produces different valid solutions?


June 20, 2013 at 01:36 am - Permalink

I was going to say that your question resembles a linear programming problem, but then I noticed who the original post was from, asking about linear programming on this forum a few months ago. So I presume formulating this as an equivalent linear programming problem and passing it to an external solver for linear programming problems is not working for you. Still, this should be the first step in any of the linear programming algorithms using the interior point method. Perhaps you can find out how it is solved there?

ChrLie wrote:
Igor wrote:

Clearly you could throw this into some optimize or root finding code but I'd be curious if there are other properties of the matrix AA that one could use. In particular, if AA is not random, I'd start by looking at its SVD.


...

Is there any way to tweak MatrixLLS, such that it produces different valid solutions?


I think what you're asking for is what AG suggested: the Singular Value Decomposition. If you call MatrixSVD to operate on the transpose of AA, then the matrix V can be interpreted as a coordinate transformation of c so that the first three (assuming matrix AA has rank 3, as in your example) rows of V^t give exactly the linear combinations of elements of c needed to solve the system, and the remaining seven rows give linear combinations of elements of c that do not affect the solution. So you can pick a number, multiply one of the last 7 rows of V^t by that number and add that to the solution to produce a different valid solution.

June 20, 2013 at 06:20 am - Permalink

Yes, the linear programming issue and this one are related, both aiming for Gibbs Free Energy minimisation. LP doesn't need a guess but calling a unix executable and first writing the problem to disk makes it very slow when calculating equilibrium e.g. along a pressure-temperature trajectory. Now I'm trying a "classical" algorithm after White, Johnson, Dantzig; J. Chem. Phys. 28, 751; 1958 which uses Lagrange multipliers. This seems to work fine and is of course must faster than the "pseudo" XOP, but it requires an initial guess.

Thanks for pointing me in the right direction. My matrix math is a bit rusty!

June 21, 2013 at 01:17 am - Permalink

Another way to approach this problem of non-negative concentrations, is to work with the logarithms of the concentrations for all the numerics. This will make things that were linear non-linear, which could be bad for the numerics, but it might make other parts of the problem more linear, and it will automatically take care of the constraints. It also avoids any numerical problems that occur if the solution is much closer to the constraint boundary than the initial guess was.

June 21, 2013 at 01:13 pm - Permalink