How can I get the derivative of a complex function including reshape, max operation?

I have a function name I want to optimize with to minimize this function

And I want to provide the gradient of function to speed up the optimization.

The definition of is

function [F] = f(x)
    % size(x) = [a*b, 1]
    Z = reshape(x, [a, b]); % size(Z）= [a, b]
    % size(A) = [c, a]
    B = A * Z; % size(B) = [c, b]
    % size(C) = [c, b]
    D = sum(B.*C, 2); % [c, 1]
    E = abs(D).^2; % [c, 1]
    F = max(E); % scalar, my objective
end

I want to use chain rules, but I have the following questions:

1. I cannot figure an easy way to get the jacobian of , i.e. the reshape function.
2. the jacobian of
3. the jacobian of function.

Thank you for your help!!

Answer 1

[Matt J]

The Jacobian of B(x) is

J1=kron(speye(b),A);

The Jacobian of D(B) is,

J2=kron(ones(1,b),speye(c));
J2(logical(J2))=C;

The max(E) function is not differentiable everywhere, so it is not clear how well fmincon will handle this problem, but at the points where it is differentiable, the Jacobian is,

 N=length(E);
 e=zeros(1,N);
 
[~,idx]=max(E);
 e(idx)=1;
 
J3=spdiags(e(:),0,N,N);

[Matt J]

As a general, but brute force tool for computing the Jacobian of linear operations, you can also use my func2mat utility, e.g.,

Jacobian = func2mat( @(Z) sum(  A*Z.*C, 2)  , ones(a,b) );

[goose neuralnet]

Thank you very much!

But sorry, I don't know how to verify this answer is true.

Another question:

1. Why the derivative of has no A in the derivative?

[goose neuralnet]

And from the definition of jacobian,

size(J1) = [a*b, b*b]
size(x) = [a*b, 1] % seems the size(J1) should be [b*b, a*b] and I can perform J1*x,
% but that's ok, I can just use x' as x.

And In my mind the derivative of a matrix w.r.t to a vector should be a tensor with shape . But I'm confused.

[goose neuralnet]

I'm interested in the theory behind your answer. Can you provide some link or matrials for me to study? Thanks!

[Matt J]

Sorry, I had a typo. The Jacobian of B(x) is

J1=kron(speye(b),A);

This is due to,

https://en.wikipedia.org/wiki/Kronecker_product#Matrix_equations

[goose neuralnet]

I seems understand the answer.

You implictly do the following things.

Faltten B into a column vector and everything seems to make sense. That's how I verify , but I still think you make mistake in .

J1 = kron(eye(b), A); % not kron(eye(b), Z)

[goose neuralnet]

Thank you!

I think I fully understand your answer!!!

But I have the following questions:

1. How can I derive these Jacobian on my own. I can verify your answer for other jacobians and I derive by myself. But I don't have a systematic way of dong this. Can you provide some material for me?
2. Although my objective F and variable x is real, but the middle consequences can be complex. So if I have complex numbers, how can I compute the jacobian?

Thank you very much!!! I'm already statisfied.

[Matt J]

Thank you very much!!! I'm already statisfied.

I'm glad, but I think my other answer is better. As you'll see there, it reaches a much more compact and efficient expression for the final gradient.

can verify your answer for other jacobians and I derive J1 by myself. Can you provide some material for me?

I already did in my earlier comments. It's just a basic Kronecker product theorem that I gave you a wiki link for.

Although my objective F and variable x is real, but the middle consequences can be complex. So if I have complex numbers, how can I compute the jacobian?

The caclulus will depend on the expression.

Answer 2

[Matt J]

Another way to see the calculation is that if

[F,i]=max(E)

then, assuming we are at a point of differentiability, F is given locally by F=D(i)^2, or,

ai=A(i,:).';
ci=C(i,:).';

F=( ai.'*Z*ci )^2

Since this is satisfied in a local neighborhood of Z, we can readily take the gradient of the expression on the right-hand side. The result, when shaped as a size [a,b] matrix, is,

Fgradient = 2*(ai.'*Z*ci)*ai*ci.'

[goose neuralnet]

Thank you very much! I'll definitely take a look at the wiki page.

What if is complex?

I know the derivative of when D is real, but I unfortunately find that my D is complex. How can I get complex jacobian?

[Matt J]

One way is to rewrite as,

and take the gradient of each term.

[goose neuralnet]

Sorry, I come back for help again.

I use fmincon with options: 'CheckGradients',true. And matlab tells me my grad is wrong. But I have checked several times and cannot figure where the problem is.

I explain my solution in the code and thank you very much!

function [F, Jacob] = f(x)
    % size(x) = [a*b, 1]
    Z = reshape(x, [a, b]); % size(Z）= [a, b]
    % size(A) = [c, a]
    B = A * Z; % size(B) = [c, b]
    % size(C) = [c, b]
    D = sum(B.*C, 2); % [c, 1]
    E = abs(D); % [c, 1]
    [F, I] = -min(E); % scalar, my objective
    
    % J1 = kron(speye(b),A); 
    % J2=kron(ones(1,b),speye(c));
    % J2(logical(J2))=C;
    
    % Assume F = E(I) for current points
    % E(I) = sqrt(real(D(I)).^2 + imag(D(I)).^2);
    % abs_square = 1/norm(E(I));
    
    % Jacobian of E(I) w.r.t to real(D(I)) is 2*real(D(I))*abs_square
    % Jacobian of E(I) w.r.t to imag(D(I)) is 2*imag(D(I))*abs_square
    
    % D = diag(A * Z * C.')
    % D(I) = A(I, :) * Z * C(I,:).'
    
    % real(D(I)) = real(ai*Z*ci.') % my ai, ci are row vector, sorry for the confusion
    % imag(D(I)) = imag(ai*Z*ci.')
    % I assure my Z = reshape(x, [a, b]) is a real matrix, hence
    % real(D(I)) = real(ai)*Z*real(ci.') - imag(ai)*Z*imag(ci.')
    % imag(D(I)) = real(ai)*Z*imag(ci.') + imag(ai)*Z*real(ci.')
    % Hence,
    % Jacobian of real(D(I)) w.r.t to x is kron(real(ci), real(ai)) - kron(imag(ci), imag(ai));
    % Jacobian of image (D(I)) w.r.t to x is kron(real(ci), imag(ai)) + kron(imag(ci), real(ai));
    
    if nargout > 1 % gradient required
        ai=tdlmat(I,:); % size(ai) = [1, J]
        ci=A_theta(I,:); % size(ci) = [1, M]
        
        J_real = kron(real(ci), real(ai)) - kron(imag(ci), imag(ai));
        J_imag = kron(real(ci), imag(ai)) + kron(imag(ci), real(ai));
        
        abs_square = 1/norm(E(I));
        Jacob = 2*real(E(I))*abs_square*J_real + 2*imag(E(I))*abs_square*J_imag;
        Jacob = Jacob.';
    end
end