Best Practices for Eigenvalue Problems

[peterjrohrer]

I am wondering if there are any established best practices or recommendations for solving eigenvalue/eigenvector problems in an OpenMDAO model and providing high-quality derivatives?

For some context: one part of my model constructs a finite element model of a structure, and solves the eigenvalue problem to find structural modeshapes and natural frequencies. The eigenmatrix $[\text{A}]$ is based on the FEM mass and stiffness matrices ( $[\text{A}] = [\text{M}]^{-1}[\text{K}]$ ) and is real but not symmetric. In the rest of the model I use a reduced basis of the first $m$ (lowest natural frequency) eigenvectors and mass-normalize the eigenvectors ( $\Phi^T [\text{M}] \Phi = [\text{I}]$ ) so there is a bit of 'post-processing' to be done after the eigenproblem to sort and select the reduced basis. For some idea of scale, the eigenmatrix has a rank of $n=500$ , and the reduced basis I am interested in is roughly rank $m=10$.

Computing and 'post-processing' the eigenvalues is straightforward using scipy.linalg.eig, though ideally I think I would use scipy.sparse.linalg.eigs to avoid calculating the full-basis in the first place. Where I am struggling a bit is calculating and checking the derivatives. For the eigenproblem, I am implementing the analytic algorithmic differentiation from He, Jonsson, and Martins, (2022) using the compute_jacvec_product method. While I see good agreement between forward and reverse mode partial derivatives under check_partials , the finite difference check is way off and seems to have quite a bit of sensitivity to step-size, making me think my step-size choice or method is not correct. Additionally, I have really reduced the size of the problem (to around $n=50$ ) to be able to run check_partials with FD in a reasonable amount of time. In the components I have written for 'post-processing' I have used the compute_jacvec_product method with algorithmic differentiation from Giles (2008) to avoid working with very large partial derivative matrices.

Are there more robust or efficient ways to include an eigenproblem like this in OpenMDAO models? Or any general tips for implementing analytic derivatives here?

This seems like a fairly common component for many different problems, so I figured it was worth asking about here rather than StackOverflow!

[JustinSGray]

Warning: I am not an expert on computing eignenvalues and the derivatives of them --- but I did stay at a holiday in express last night :)

#######################################################

Broadly, there are two approaches I have seen for dealing with egienvalues in a differentiable context.

1. Compute them, then pass them all through some kind of aggregation function
2. Compute them, and keep them separate --- but add some kind of methods to avoid mode switching or otherwise manage it via mode tracking

In either case, you still need a method to compute eigen values and their derivatives. In my (just enough knoweldge to be dangerous) opinion, the method that He, Jonsson, and Martins offers is one of the better approaches for that in an aero-structural context. Graeme Kennedy also has some numerical approaches that he's developed for frequency constraints in topology optimization as well, which are integrated into newer versions of TACS --- and thus might be able to be exposed to OpenMDAO in a more direct way.

In both cases, I am pretty confident that scipy's eigenvalu solver wasn't up to the task of the size of problems that He et. al. and Kennedy were solving. If your problem is small enough to be within the scope of those algorithms though, it understandable why you might want to go that way.

#########################################################

Its not even a little surprising to me that your FD derivatives are way off. The eigenvalue approximation methods are noisy, and its really tough to make them small-scale-stable for very stiff matricies like you would see in an FEM. I have worked on other problems where I simply could not get accurate/stable finite difference derivatives at all (mainly pyCycle problems). So in this situation, if you can't use something like CS as an independent verification ... i might let the optimization itself indicate if the derivatives are "correct" (maybe I should say "correct enough"). If the optimization works, then you're good.

[markleader]

I'll start by saying I agree with everything Justin suggested, and emphasize using aggregation on your eigenvalues to prevent derivative discontinuities due to mode-switching.

Additionally, here are a couple of techniques that I found useful for natural frequency constraints with topology optimization:

1. Using a Jacobi-Davidson solver as opposed to a Lanczos solver; the Jacobi-Davidson method solves to a looser tolerance (say O(10^-6)) than Lanczos (which might be O(10^-12)), but for optimization this is totally fine and saves you a ton of time.
2. Using eigenvalue recycling; use the converged eigenvector from the previous optimization iteration as the starting point in the eigenvalue solve for the current optimization iteration. Since the structure is generally not changing drastically between optimization iterations, this gives you a really close initial estimate and improves the speed of the eigenvalue solve significantly.

These are probably only necessary for really large problems, but worth noting as far as speed improvements go. This paper (Section 3B) has more of the details if you are interested.

If you try Justin's suggestions and are still having trouble, I would suggest trying TACS as your finite element solver, which has an efficient and robust modal solver implemented with derivatives available as well. You can put TACS into an OpenMDAO component and connect it to your existing OpenMDAO problem. I would be happy to help if you are interested in this route.

[peterjrohrer]

Thank you both for the replies! Though I'm not sure I follow what you mean by applying an aggregation function or mode tracking in this context?

I think we will look into using TACS for now to avoid re-inventing the wheel, but I am of course also open to other suggestions.