How do I minimize a particular error function with lsqnonlin?

Hello all,

I am trying to minimize the error function

E = abs(vecnorm(b_hat) - vecnorm(b_ref))

by solving for A and x in

b_hat = A * (b_raw - d)

where b_ref, b_raw and b_hat have dimensions of 3 x 100 (they are all 3 dimensional vectors, with 100 different values over time). A is a 3 x 3 matrix and d is a 3 x 1 vector. Only the values of b_raw and vecnorm(b_ref) are known to me. The correct values of A and d would help me to calculate b_hat, whose norm would be equal to that of b_ref. At the moment, I am solving them using lsqnonlin with the option 'trust-region-reflective', which helps me also to set up the upper and lower bounds, like this:

 fun = @(x) abs(vecnorm([x(1) x(2) x(3); x(4) x(5) x(6); x(7) x(8) x(9)] * (b_raw - [x(10); x(11); x(12)])) - vecnorm(b_ref));
        x0(1, 1:9) = [1 0 0 0 1 0 0 0 1];
        x0(1, 10:12) = [0  0  0];
        
        options = optimoptions('lsqnonlin');
        
        options.Algorithm = 'trust-region-reflective';
        lb = [0.8 -0.1 -0.1 -0.1 0.8 -0.1 -0.1 -0.1 0.8 -100000 -100000 -100000];
        ub = [1.2 0.1 0.1 0.1 1.2 0.1 0.1 0.1 1.2 100000 100000 100000];
        
        options.StepTolerance = 1e-10;
        options.FunctionTolerance = 1e-10;
        options.OptimalityTolerance = 1e-10;
        options.MaxFunctionEvaluations = 1e10;
        options.MaxIterations = 1e3;
        
        [y, resnorm] = lsqnonlin(fun, x0, lb, ub, options);
        A = [y(1) y(2) y(3);...
             y(4) y(5) y(6);...
             y(7) y(8) y(9)];
        d = [y(10); y(11); y(12)];

Is there a better way to solve this? If I set the limits far-apart too much, the estimation of A and d becomes incorrect.

Thanks a lot!

[Matt J]

So, you are deliberately looking for a weaker condition in your solution than A*b_raw-d ≈ b_ref?

Also, is A supposed tobe a rotation matrix?

[Debdeeprc]

Hi Matt,

Thanks for your comment. Please note, A * b_raw-d is not even close to b_ref. Only, their norms are "almost equal". I will give more context here to clarify. b_raw is actually the magnetic field output vector from an uncalibrated magnetic field sensor which is on a low earth orbiting satellite. The sensor output values (b_raw) have to be corrected using a 3 x 3 calibration matrix A and an offset vector d to get the calibrated output b_hat. The information about what the correct magnetic field norm should be is already available from the IGRF magnetic field model, which in my case depends on the satellite's position. The IGRF gives the output in an Earth Centred frame (b_ref). A calibrated magnetic field sensor output should have the same norm/magnitude as b_ref, however the individual vector elements would be different as they are in two different reference frames: b_hat is in satellite body fixed frame and b_ref is Earth centred.

Thanks again!

[Debdeeprc]

Hi Matt,

A is NOT a rotation matrix. It is more commonly referred to as a calibration matrix, whose diagonal terms scale the raw sensor values and the off-diagonal terms correct the orthogonality errors in the magnetic field sensor (sometimes cheaper sensors don't have their axes mutually perpendicular to one another). Actually, you can read more about this specific part under the section Correction Technique of https://www.mathworks.com/help/nav/ug/magnetometer-calibration.html

Thanks!

Answer 1

[Matt J]

I don't see much that I think can be improved, except perhaps to formulate it so that your errors are differentiable - the lsqnonlin algorithms are generally driven by Jacobian calculations.

v_ref=vecnorm(b_ref).^2; %compute only once

B=[b_raw; -ones(1,100)]; %put into homogeneous form

fun = @(Ad) vecnorm(Ad*B).^2 - v_ref; %Use difference of vecnorm^2 for differentiability

         options = optimoptions('lsqnonlin');
        
        options.Algorithm = 'trust-region-reflective';
        lb = [0.8 -0.1 -0.1 -0.1 0.8 -0.1 -0.1 -0.1 0.8 -100000 -100000 -100000];
        ub = [1.2 0.1 0.1 0.1 1.2 0.1 0.1 0.1 1.2 100000 100000 100000];
        
        options.StepTolerance = 1e-10;
        options.FunctionTolerance = 1e-10;
        options.OptimalityTolerance = 1e-10;
        options.MaxFunctionEvaluations = 1e10;
        options.MaxIterations = 1e3;
        
        [Ad, resnorm] = lsqnonlin(fun, eye(3,4), lb, ub, options);
        
 A=Ad(1:3,1:3);
 d=Ad(:,4);

[Debdeeprc]

Hi Matt,

Thanks a lot for your suggestions. I already see improvements in the results after implementing your suggestions, especially in cases where sufficient variation in vecnorm(b_ref) is not available. Using the difference of vecnorm^2 instead of vecnorm works better.

Have a great day,

Debdeep