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

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:

[Torsten]

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]

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...

Answer 1

[Matt J]

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

1 0 0

I have submitted a bug report.

[Matt J]

The bad logic appears to be in this section:

if ishermitian(A) || ishermitian(A, 'skew')
    [Q, d] = eig(full(A), 'vector');
    expd = exp(d);
    F = (Q.*expd.')*Q';
    if isreal(expd)
        F = (F+F')/2;
    end
    return;
end

I think it should really be,

if ishermitian(A) || ishermitian(A, 'skew')
    [Q, d] = eig(full(A), 'vector');
    expd = exp(d);

    F = real( (Q.*expd.')*Q'  );
    
    return;
end

[Torsten]

What if A is complex ?

A = [1 3-1i 4;3+1i -2 -6+1i;4 -6-1i 5];
ishermitian(A)

ans = logical
   1

[Q, d] = eig(full(A), 'vector');
d

d = 3×1

-7.9622 2.7363 9.2259

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

expd = exp(d);
F = real( (Q.*expd.')*Q'  );
F

F = 3×3

1.0e+03 * 0.9612 -1.0241 2.5403 -1.0241 1.7128 -3.4018 2.5403 -3.4018 7.4982

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

●

expm(A)

ans = 
   1.0e+03 *

   0.9612 + 0.0000i  -1.0241 - 0.7486i   2.5403 + 0.8175i
  -1.0241 + 0.7486i   1.7128 + 0.0000i  -3.4018 + 1.1120i
   2.5403 - 0.8175i  -3.4018 - 1.1120i   7.4982 + 0.0000i

[Matt J]

True. So, I guess it should be...

if ishermitian(A) || ishermitian(A, 'skew')
    [Q, d] = eig(full(A), 'vector');
    expd = exp(d);
    F = (Q.*expd.')*Q';
    if isreal(A)
        F = real(F);
    end
    return;
end

[Matt J]

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 =

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: []

[xingxingcui]

@Matt J Thank you for your keen insight!

[Paul]

@Matt J

Great catch. I really hope that error gets fixed before the 2026A is formally released into the wild. Nevertheless, seems strange and, frankly, quite disappointing that it made it even into a pre-release.

Can you update your comment above to use format long to show expm(A) to full precision?

[Matt J]

@Paul I updated it, but even format long isn't enough to reveal the non-zeroness of the imaginary part.

[Paul]

Now I'm confused. If I run the bad logic here I get the expected (real) result, but the same code in 2026a yields a complex result?

A=zeros(3); A(3,2)=1; A(2,3)=-1;
[Q, d] = eig(full(A), 'vector');
expd = exp(d);
F = (Q.*expd.')*Q'

F = 3×3

1.0000 0 0 0 0.5403 -0.8415 0 0.8415 0.5403

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

expm(A)

ans = 3×3

1.0000 0 0 0 0.5403 -0.8415 0 0.8415 0.5403

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Matt J]

If I run the bad logic here I get the expected (real) result, but the same code in 2026a yields a complex result?

On my computer, in both R2024b and R2026a,:

F =

1.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i

0.0000 + 0.0000i 0.5403 - 0.0000i -0.8415 - 0.0000i

0.0000 + 0.0000i 0.8415 + 0.0000i 0.5403 - 0.0000i

However, I do not have R2025a installed, so I cannot compare, but a real-valued result is theoretically unstable numerically, when the eigenvalues are complex, as they are here. The differences we see would have to be floating point noise -- possibly a CPU-dpendent thing, or version-dependent differences in multithreading.

[Paul]

@Matt J

I'm running 2024a and your A matrix causes expm to go down a different path because in 2024a expm checks ishermitian(A) but does not check for ishermitian(A,'skew').

In 2024a we have:

if ishermitian(A)
    [Q, d] = eig(full(A), 'vector');
    F = (Q.*exp(d).')*Q';
    F = (F+F')/2;
    return;
end

from which I gather that if A is Hermitian then F must be conjugate symmetric. But I don't think your latest proposal here would ensure that.

I assume that the skew check and other mods were added in 2024b, which is why you're seeing the complex result? I don't have 2024b handy so can't check myself.

But @郦歌 reports no issues in 2025b and prior, albeit for different input matrices than we're playing with here.

[Matt J]

F = (Q.*exp(d).')*Q' is the exact formula for exp(A) so any failure of F to have requisite properties like realness or Hermitian-ness would have to be due to floating point noise. You can make sure the floating point representation of the final result is ideally Hermitian by modifying as follows.

isherm=ishermitian(A);

if isherm || ishermitian(A, 'skew')
    [Q, d] = eig(full(A), 'vector');
    expd = exp(d);
    F = (Q.*expd.')*Q';
    
    if isreal(A)
        F = real(F);
    end
    
    if isherm
        F=(F+F')/2;
    end
    
    return;
end

[Paul]

According to Wikipedia, the properties that hold are:

"If X is symmetric then eX is also symmetric, and if X is skew-symmetric then eX is orthogonal. If X is Hermitian then eX is also Hermitian, and if X is skew-Hermitian then eX is unitary."

I don't think the else case above is correct. And the real(F) should be predicated on A being real.

For the cases where F is supposed to be orthogonal or unitary, what would be the best way to ensure that?

[Matt J]

I think that's what I now have.

[Matt J]

For the cases where F is supposed to be orthogonal or unitary, what would be the best way to ensure that?

I don't think you can. Floating point Hermitian-ness is something that can be ensured (as above), but probably not orthogonality. F*F' = I will always give some floating point error.

[Paul]

One more sublety that I just realized. The isreal check will fail if the imaginary part of a complex matrix zero.

A = [0,0,0;0,0,-1;0,1,0];
A = complex(A,0*A);
ishermitian(A,'skew')

ans = logical
   1

isreal(A)

ans = logical
   0

Perhaps need to change to

if ~any(imag(A),'all')
    F = real(F);
end

Unrecognized function or variable 'F'.

[dpb]

"The isreal check will fail if the imaginary part of a complex matrix zero."

Yes, isreal only checks the data class of the variable; it doesn't care what the values are.

And then you're back into the problem of how much is sufficient that an array is intended to be real but has floating point rounding in the imaginary part making it complex? Just as the earlier comment that (almost) always the ==1 will fail for checking unity. Floating point is a b@!xx! when it comes to such details. Numerical computing and math rarely overlap neatly.

[Matt J]

@dpb I think @Paul was suggesting that if A has an imaginary part, whether explicit or implicit, that is perfectly zero then theoretically so should the output of exp(A). So, you might do something like this,

isherm=ishermitian(A);

if isherm || ishermitian(A, 'skew')
    
    [Q, d] = eig(full(A), 'vector');

    F = (Q.*exp(d).')*Q';
    
    if isreal(A)
      F = real(F);
    elseif ~any(imag(A(:)))  %real, but of complex type
      F=complex(real(F), imag(A));
    end
    
    if isherm
        F=(F+F')/2;
    end
    
    return;
end

[dpb]

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]

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