# Wrong transpose leads to complex step gotcha, matlab

## weblog/

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: - it's real-valued on the real line ✅ - evaluated at points in the complex plane near the input ✅ - holomorphic

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