When translating math to matlab, I admit I usually translate xᵀ into x'. For example, suppose x∈ℝ⁴ and I have a quadratic energy:
xᵀ A x
I would implement this in matlab with
x' * A * x
The ' isn't actually calling transpose. It's actually calling ctranspose, the complex conjugate transpose. For real input, this is the same as transpose so most of the time this works just fine for me.
However, today I was experimenting with the complex step method for numerically computing gradients. I immediately got incorrect results. Here's a simple example:
i = sqrt(-1);
t = 1.5;
A = [1 2 3;2 3 4;3 4 5];
x = @(t) [1;t;t^2];
f = @(t) x(t)' * A * x(t);
fd = (f(t+1e-5)-f(t-1e-5))/(2*1e-5);
cs = imag(f(t + i*1e-5))/1e-5;
fprintf('finite difference: %g\n',fd);
fprintf('complex step: %g\n',cs);
This prints:
finite difference: 152.5
complex step: 0
The trouble is the x(t)'. My f function is supposed to be acting only on the real inputs, so that the imaginary parts can "track" the first derivative. Unfortunately, the ctranspose mangles this and mixes in its derivatives with respect to the imaginary parts.
As wikipedia writes, the requirement on f to use complex step method is that:
The failure to be holomorphic might come as a surprise if you casually think of it as "differentiable and innocuous looking". Certainly my function looks simple enough. But because my f uses complex conjugation, it became non-holomorphic. We can see that complex conjugation is not holomorphic by looking at the Cauchy-Reimann equations:
Suppose f(x + yi) = u(x,y) + v(x,y)i = x - yi. Then:
∂u/∂x = 1 ≠ ∂v/∂y = -1.
The fix is simple. In Matlab the transpose is written .'. So I just change x(t)' above to x(t).' and it prints
finite difference: 152.5
complex step: 152.5
note: as far as I can tell, using standard finite differences to double check if a function satisfies the Cauchy-Riemann equations near t basically amounts to just checking if complex-step agrees with standard finite differencing:
x = t;
y = 0;
dudx = real( f(x+1e-5 + y*i) - f(x-1e-5 + y*i) ) / (2*1e-5);
dvdy = real( f(x + (y+1e-5)*i) - f(x + (y-1e-5)*i) ) / (2*1e-5);
abs(dudx - dvdy)
So, if I'm worried my function f uses wrong transposes I guess I either need to search for those or resort to checking if standard finite differencing matches complex step for the inputs I care about. Is there a more bullet proof way to check?
Update: Maybe something like this could be useful. This is a class that extends double which throws an error if you call ' on it. Save it in a file NoCTranspose.m:
classdef NoCTranspose < double
methods
function obj = ctranspose(data)
error("You called ctranspose (') on a NoCTranspose object");
end
end
end
Then if you try to use it using the original code above:
t = NoCTranspose(1.5);
f(t)
you get an error:
Error using '
You called ctranspose (') on a NoCTranspose object
Error in test_complex_step>@(t)x(t)'*A*x(t) (line 29)
f = @(t) x(t)' * A * x(t);
update: another complex step gotcha in matlab, min(X,Y) will not work correctly because if X is complex it defaults to using abs (a reasonable choice but we want to branch based on the real value) so you have to use min(X,Y,‘ComparisonMethod’,‘real’);