Why the optimization results of lsqnonlin are different in R2026a and R2025a?

187 次查看（过去 30 天）

郦歌 2025-11-10，14:21

2
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2181154-why-the-optimization-results-of-lsqnonlin-are-different-in-r2026a-and-r2025a

评论： Matt J 2025-11-13，22:35

I used the following code to solve a nonlinear optimization problem. All random seeds and optimization options are the same, but the result in the R2026a pre-release is complex-valued, whereas in R2025a or older releases the result is real-valued. All code and data are attached.

% R2025b or earier is real-base solution
% R2026a is complex-base solution?

R2025b:

load temp.mat
rng default;
options = optimoptions('lsqnonlin', 'Algorithm','levenberg-marquardt', 'Display','final',...
    'MaxFunEvals',2000, 'MaxIter',1e3, 'TolFun',1e-6, 'TolX',1e-6, 'Jacobian','off');
[x,resnorm,~,exitflag,output] = lsqnonlin(@(p)residual_KR_robust(double(X1_ok), double(X2_ok), imsize1, imsize2, p, 1000), [1000 1000 0 0 0],...
    [],[])
Solver stopped prematurely.

lsqnonlin stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 5.000000e+02.
x = 1×5
1.0e+04 *

    2.3371    2.4754    0.0000   -0.0000    0.0000
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
resnorm = 2.8483e+07
exitflag = 0
output = struct with fields:
      firstorderopt: 1.6500e+07
         iterations: 83
          funcCount: 504
       cgiterations: 0
          algorithm: 'trust-region-reflective'
           stepsize: 1.2800e+03
            message: 'Solver stopped prematurely.↵↵lsqnonlin stopped because it exceeded the function evaluation limit,↵options.MaxFunctionEvaluations = 5.000000e+02.'
       bestfeasible: []
    constrviolation: []

R2026a:

2 个评论
显示无隐藏无

Torsten 2025-11-10，15:44

Starting from real-valued initial guesses for the parameters, do you have a line in your code that could produce complex results for an element of "err" ? I cannot find one - thus I have no explanation why you could end up with complex-valued parameters at all.

dpb 2025-11-10，16:14

编辑：dpb 2025-11-10，16:36

在 MATLAB Online 中打开

type residual_KR_robust
function [ err ] = residual_KR_robust( X1, X2, imsize1, imsize2, paras, sigma )

% parameretes sigma indicates the distance scope for inliers with homopraphy matrix, in pixels

k1 = paras(1);
k2 = paras(2);
theta = paras(3:5);
% yaw = paras(3);
% pitch = paras(4);
% roll = paras(5);

K1 = [k1, 0, imsize1(2)/2;
     0, k1, imsize1(1)/2;
     0,  0, 1];
K2 = [k2, 0, imsize2(2)/2;
      0, k2, imsize2(1)/2;
      0,  0, 1];
theta_m = [0         -theta(3) theta(2)
           theta(3)  0         -theta(1)
           -theta(2) theta(1)  0];
R = expm(theta_m);
% Ry = [cos(yaw),     0,              -sin(yaw);
%       0,            1,              0;
%       sin(yaw),     0,              cos(yaw)] ;
% Rp = [1,            0               0;
%       0,            cos(pitch),     sin(pitch);
%       0,            -sin(pitch),    cos(pitch)] ;
% Rr = [cos(roll),    sin(roll),      0;
%       -sin(roll),   cos(roll),      0;
%       0,            0,              1] ;
% R = Rr * Rp * Ry;

err = residual_H(X1, X2, K1, K2, R);

outlier = (abs(err) > sigma);
err(outlier) = sign(err(outlier)) .* (sigma + sigma * log(abs(err(outlier))/sigma));
% err(outlier) = sign(err(outlier)) .* sqrt(2*sigma*abs(err(outlier)) - sigma*sigma);

end
type residual_H
function [ errH ] = residual_H( X1, X2, K1, K2, R)

% homopraphy matrix
H = K2*R/K1;

X1_p2 = H * X1;
X1_p2(1,:) = X1_p2(1,:) ./ X1_p2(3,:) ;
X1_p2(2,:) = X1_p2(2,:) ./ X1_p2(3,:) ;
errH1_2 = X2 - X1_p2;
X2_p1 = H \ X2;
X2_p1(1,:) = X2_p1(1,:) ./ X2_p1(3,:) ;
X2_p1(2,:) = X2_p1(2,:) ./ X2_p1(3,:) ;
errH2_1 = X1 - X2_p1;

errH = [errH1_2(1,:)', errH1_2(2,:)', errH2_1(1,:)', errH2_1(2,:)'];

end

@Torsten asked "do you have a line in your code that could produce complex results...?"

The backslash solve might under some circumstances, couldn't it?

Setting a breakpoint on condition ~isreal(errH) would be able to discover when it first occurs and let figure out what might have happened and (maybe then) why.

I don't suppose there's anything about lsqnonlin in the release notes...

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Matt J 2025-11-10，20:31

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2181154-why-the-optimization-results-of-lsqnonlin-are-different-in-r2026a-and-r2025a#answer_1571810

编辑：Matt J 2025-11-11，2:53

There appears to be a new (and buggy) implementation of expm.m in R2026a, resulting in incorrectly complex results for skew symmetric matrices:

>> A=zeros(3); A(3,2)=1; A(2,3)=-1

A =

0 0 0

0 0 -1

0 1 0

>> B=expm(A)

B =

1.000000000000000 + 0.000000000000000i 0.000000000000000 + 0.000000000000000i 0.000000000000000 + 0.000000000000000i

0.000000000000000 + 0.000000000000000i 0.540302305868140 - 0.000000000000000i -0.841470984807896 - 0.000000000000000i

0.000000000000000 + 0.000000000000000i 0.841470984807896 + 0.000000000000000i 0.540302305868140 - 0.000000000000000i

>> imag(B)==0

ans =

3×3 logical array

1 1 1

1 0 0

I have submitted a bug report.

19 个评论
显示 17更早的评论隐藏 17更早的评论

Matt J 2025-11-11，0:33

编辑：Matt J 2025-11-11，0:38

在 MATLAB Online 中打开

This fixes the lsqnonlin result. Of course, I don't know how much you're inclined to trust the rest of the pre-release with the built-in expm out of whack...

function [ err ] = residual_KR_robust( X1, X2, imsize1, imsize2, paras, sigma )
k1 = paras(1);
k2 = paras(2);
theta = paras(3:5);
ctr1=imsize1/2;  %NOTE: These are the wrong expressions for the center of a 1-based
ctr2=imsize2/2;  %      image coordinate grid. Should be ctr = (imsize+1)/2
K1 = [k1, 0, ctr1(2);
     0, k1,  ctr1(1);
     0,  0, 1];
K2 = [k2, 0, ctr2(2);
      0, k2, ctr2(1);
      0,  0, 1];
theta_m = [0         -theta(3) theta(2)
           theta(3)  0         -theta(1)
           -theta(2) theta(1)  0];
[V,D]=eig(theta_m);             %<-----Matt J
R=real(V*diag(exp(diag(D)))/V); %<-----Matt J
err = residual_H(X1, X2, K1, K2, R);
outlier = (abs(err) > sigma);
err(outlier) = sign(err(outlier)) .* (sigma + sigma * log( abs(err(outlier))  / sigma )    )  ;
end

Solver stopped prematurely.

lsqnonlin stopped because it exceeded the function evaluation limit,

options.MaxFunctionEvaluations = 5.000000e+02.

x =

1.0e+04 *

2.3371 2.4754 0.0000 -0.0000 0.0000

resnorm =

2.8483e+07

exitflag =

output =

struct with fields:

firstorderopt: 1.6500e+07

iterations: 83

funcCount: 504

cgiterations: 0

algorithm: 'trust-region-reflective'

stepsize: 1.2800e+03

message: 'Solver stopped prematurely.↵↵lsqnonlin stopped because it exceeded the function evaluation limit,↵options.MaxFunctionEvaluations = 5.000000e+02.'

bestfeasible: []

constrviolation: []

dpb 2025-11-11，18:55

在 MATLAB Online 中打开

I was just commenting back on @Paul wrote as

if ~any(imag(A),'all')

will only be true if all(imag(A))==0 identically; even the LSB in one element will fail such that if A got converted from intended real to complex owing to some computational gaff like under discussion here, should it really be complex with rounding error complex part or real?

Matt J 2025-11-13，22:35

Tech support has informed me that they are aware of the issue, and that it will be fixed in the R2026a general release.

请先登录，再进行评论。

类别

Mathematics and Optimization Symbolic Math Toolbox MuPAD Mathematics Equation Solving Numeric Solvers

在 Help Center 和 File Exchange 中查找有关 Numeric Solvers 的更多信息

产品

Optimization Toolbox

Community Treasure Hunt

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

Start Hunting!

Translated by

Why the optimization results of lsqnonlin are different in R2026a and R2025a?

2 个评论
显示无隐藏无

采纳的回答

19 个评论
显示 17更早的评论隐藏 17更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

Community Treasure Hunt

Why the optimization results of lsqnonlin are different in R2026a and R2025a?

2 个评论 显示 无隐藏 无

采纳的回答

19 个评论 显示 17更早的评论隐藏 17更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

Community Treasure Hunt

2 个评论
显示无隐藏无

19 个评论
显示 17更早的评论隐藏 17更早的评论