Parameter covariance from lsqnonlin when using Jacobian Multiply function?

Hi, I'm trying to estimate the uncertainty in my model parameters for a large NLSQ problem (~2M measurements, ~10k unknowns), so I'm constantly struggling to manage RAM even with 128GB on my machine. I have the algorithm working reasonably well using the Jacobian Multiply function interface, but now when it returns I'd like to estimate the parameter covariance in order to assess the fit and the uncertainties. Normally, if one does:

[xCurrent,Resnorm,FVAL,EXITFLAG,OUTPUT,LAMBDA,JACOB] = lsqnonlin(problem);
xCov = inv(JACOB.'*JACOB)*Resnorm/(numel(FVAL)-numel(xCurrent));

You can get the covariance, or you can just send the Jacobian into nlparci() and get confidence intervals. However, if I try that, I get "out of memory" errors, which is why I'm using the JM functions in the first place. I tried sending in eye(numel(xCurrent)) with flag zero to get the JM function to compute the inner product of the Jacobian, but same out of memory error. Is there a loop method I should use to build up the rows/columns of J.'*J? Do I try using a sparse identity matrix? I notice the JACOB returned is a sparse matrix class, but I'm not sure how sparse my native Jacobian really is (some cases moderately sparse, other times it can be pretty full rank). Is there a lower-footprint means of building up the covariance?

[Martin Ryba]

Following up on myself, a brute force way to run the JacobMult function to get individual columns of (J'J) looks like this:

nParms = numel(xCurrent);
nPts = numel(FVAL);
xCov = zeros(nParms);
for icol = 1:nParms
    Yin = zeros(nParms,1);
    Yin(icol) = 1;
    xCov(:,icol) = JMfun(xCurrent,Yin,0);
end
xCov = inv(xCov)*Resnorm/(nPts-nParms);

Sending an input vector that has only one nonzero element likely wastes a lot of operations, so the better alternative would be to build a subfunction that emulates JMpos with the single nonzero input and sends that result straight into JMneg.

[Torsten]

And you are sure it makes sense to estimate ~10k model parameters, even with ~2M measurement data ? And how to evaluate xCov in the end ? And how to generate senseful initial values for the parameters ? I'm sceptical ...

[Bruno Luong]

@Torsten this kind of dimensions is not unusual for climate/meteorogical models. I worked with similar dimension problem more than twenty year ago already.

[Torsten]

And do parameters in such models have physical meaning or are they just used to give a good fit - value and meaning unimportant ?

[Bruno Luong]

In my case of study the parameters are the initial state condition of the climate model and the data are satelite observation data , so yes the parameters have very much physical meaning.

[Guillermo Soriano]

Goog evening

I worked with an algorithm for solving Non Linear Least Square and without any problem.. but recently appears for this staement : inv(J'*inv(R)*J)*J'*inv(R)*INN to better use \ which is faster. I did and worked fine until yesterday...

today appears this message : Error using \

Matrix dimensions must agree.

The matrix are:

Name Size Bytes Class Attributes

J 450x4 14400 double

R 450x450 1620000 double

Does someone explain to me which is the right notation?

Thank you in advance

[Bruno Luong]

@Guillermo Soriano please post your code using "\".

otherwise try

lscov(J, INN, R)

Example with toy data et 4 ways to computed that quantity:

m = 10;
m = 10;
n = 4;
J=rand(m,n);
W=randn(m);
R=W*W';
INN=rand(m,1);

x1 = inv(J'*inv(R)*J)*J'*inv(R)*INN

x1 = 4×1

-1.3925 0.7634 1.0429 0.0057

K = J'/R; x2 = (K*J)\(K*INN)

x2 = 4×1

-1.3925 0.7634 1.0429 0.0057

U = chol(R); x3 = (U'\J)\(U'\INN)

x3 = 4×1

-1.3925 0.7634 1.0429 0.0057

x4 = lscov(J, INN, R)

x4 = 4×1

-1.3925 0.7634 1.0429 0.0057

[Guillermo Soriano]

Bruno Luong: hi.. thank you so much. I tryed with lscov and works very nice.. also I tryed with K = J'/R and also works fine.

Find atached part of my algorithm that process the non lineal least square for target parameters and the use of Newton Raphson to find the root. I need to reduce so many lines of code.. could be possible to use function lsnlin ? I was trying to use but without results due mainly the covariance.

Once more thnak you so much

[Bruno Luong]

I don't kwow lsnlin, but lsqnonlin.

As the method for x3 suggested, doing the model to fit data y with with covariance R

f(x) = y

is equvlalent to do the stadard least square on

U'\f(x) = U'\y

where U is the cholesky (or square root) of R.

In other word, use lsqnonlin to solve g(x) = z with

g(x) := U'\f(x)
z := U'\y?

Answer 1

[Torsten]

Can you form the (sparse) matrix

B := [JACOB -speye(size(JACOB,1)) ; zeros(size(JACOB,2)) JACOB.']

?

Then you could try to solve

JAC.'*y = e_i
JAC*x = y_i

with e_i as i_th column of eye(numel(xCurrent)) using an iterative method.

If you concatenate the y_i to a matrix, you should get the inverse of JACOB.'*JACOB.

[Martin Ryba]

No, turns out when you use the JM functions, the JACOB returned by lsqnonlin is useless, it's the placeholder sparse matrix you create in the primary fit function. See my update above regarding constructing the real J.

Answer 2

[Matt J]

Some people will compute just the diagonal D of J.'*J, approximating it as a diagonal matrix. This calculation is simply,

D=sum(J.^2)

or

D=(J.^2).'*ones(N,1)

Since you have a routine to compute J.'*x, maybe it would be similarly easy for you to calculate (J.^2).'*x;

[Bruno Luong]

But the diagonal of inverse if J'*J is more of interest.

[Matt J]

Well, yes, but if J'*J is approximately diagonal, then diag(inv(J'*J)) is approximately D.^-1

You can possibly improve the accuracy of the inversion by computing some of the off-diagonals, e.g.,

D_k = sum( J(:,1:end-k).*J(:,k:end) ); %equivalent to diag(J.'*J,k)

[Bruno Luong]

Really?

H=magic(6)

H = 6×6

35 1 6 26 19 24 3 32 7 21 23 25 31 9 2 22 27 20 8 28 33 17 10 15 30 5 34 12 14 16 4 36 29 13 18 11

D=diag(diag(H))

D = 6×6

35 0 0 0 0 0 0 32 0 0 0 0 0 0 2 0 0 0 0 0 0 17 0 0 0 0 0 0 14 0 0 0 0 0 0 11

diag(inv(H))

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 8.001607e-19.

ans = 6×1

1.0e+15 * 2.2518 -0.0000 1.1259 2.2518 0.0000 1.1259

diag(inv(D))

ans = 6×1

0.0286 0.0312 0.5000 0.0588 0.0714 0.0909

[Matt J]

I don't see the analogy. The assumption was that H=J.'*J is approximately diagonal, which doesn't appear to be the case in this example.

[Bruno Luong]

The point is what you claim inv(D) ~ inv(H) has no foundation whatsoever.

[Matt J]

It is well-founded if D=H and approximately founded if .

Answer 3

[Bruno Luong]

A low rank approximation of

H=inv(J'*J);

is using svds

r = 10;
[~,S,V]=svds(J,r,'smallest');
r = size(S,1); % == 10
s = diag(S);
G = V*spdiags(1./s.^2,0,r,r)*V';

In order to have G ~ H, you need to set r so that the partial sum of (1/s.^2), where s is the singular values of J, is close to total sum. So depending on the spectrum of your Jacobian it can be hard or not to well approximate with low value of r.

But probably the singular vector V(:,j) give plenty of information about the sensitivity of the parameters.