43 views (last 30 days)

Show older comments

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
on 23 Feb 2021

Edited: Matt J
on 24 Feb 2021

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

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!