How to calculate the jacobian of a linear symbolic system for optimization

Hi,

I have a very simple linear system:

• T’ = A*T + B*U
• Y = C*T + D*U

I try to minimize Y (which is simply equal to T(1)) with lsqnonlin solver, so in each time step the objectif function (function whose sum of squares is minimized) calculate Y by using a new diagonal base:

• [P,F] = eig(A)
• X’ = F*X + E*U
• T = P*X

with F is the diagonal matrix in the new base.

In order to help the solver, I would like to provide it with the jacobian, so the objectif function will return a second output argument with the Jacobian value.

The function Jacobian in Matlab allows me to calculate it very easily, but the concern is that the diagonalization of the system in symbolic variables does not work beyond a small size of matrix B (P\B or inv(P)*B, or linsolve (P, B)).

Does anyone have an idea how could I solve this system with symbolic tools, in order to calculate the Jacobian for the solver ?

Thank you in advance for your feedback.

Answer 1

[Walter Roberson]

"small size of matrix B (P\B" implies to me that P and B have the same number of rows. You have calculated P through eig(A) so P will have the same number of rows as A.

If your A is symbolic, to get symbolic P, then you would be calculating the symbolic eigenvectors of A.

The maximum feasible size of symbolic matrix to eig() is 4 x 4, which is not exactly fast but only takes a couple of minutes.

Taking the symbolic eig() of a 5 x 5 takes many hours. And the calculation is pretty useless as it involves roots of a general quintic, but general quintics do not have expressible roots in closed forms.

[Michael cohen]

Ok I understand, thank you for your answer,

And yes you're right : A, B and P actually have the same number of lines. And E is just equal to P \ B.

eig (A) takes indeed a long time when matrix A has more than 4 rows. For 5 rows it's still reasonable in my case, perhaps because my matrix A contains a lot of zero values (?). Anyway, the main concern comes rather during the calculation of E, and it is the calculation of the inverse of P which is extremely long.

Do you think there would be another way to calculate the Jacobian of the function or should I forget the idea of providing it to lsqnonlin?