Implementation of solver of generalized eigenvalue problem

[tky823]

Summary

In IP2, we solve the following generalized eigenvalue decomposition for positive definite matrices A and B. A @ x = lambda @ B @ x

1. Current implementation By using Cholesky decomposition L @ L.H = B, solve the following eigenvalue problem. C @ y = lambda @ y, where C = L^{-1} A @ L.H^{-1} and y = L.H @ x Then, x = L.H^{-1} @ y.

2. Implementation referred to [1] Compute G = B^{-1} @ A (not positive definite). lamb_max = (tr(G) + sqrt(tr(G)*2 - 4 det(G))) / 2 lamb_min = (tr(G) - sqrt(tr(G)*2 - 4 det(G))) / 2 x_max = (-G_12, G_11 - lamb_max) x_min = (G_22 - lamb_min, - G_21) The order in [1] may be incorrect.

Which is faster?

[1] paper of ISS2

[tky823]

experimental notebooks: https://colab.research.google.com/github/tky823/ssspy/blob/5954d98d8a3755799a24d218143b1981c81a4175/notebooks/ISSUE/issue199.ipynb

1. Current implementation

   • Pros. More stable than 2
   • Cons. Slower than 2

2. Implementation referred to [1]

   • Pros. About 15% faster than 1
   • Cons. Unstable at stationary points of cost function

[tky823]

In terms of computational stability, I will choose the current implementation.

[tky823]

I found this is due to the stability of np.linalg.eigh. The naive implementation of np.linalg.eigh for a 2x2 matrix using closed form is unstable. https://colab.research.google.com/github/tky823/ssspy/blob/bf8bd1187543ab163f8ff445040ee9e38ecb1c8a/notebooks/ISSUE/issue199.ipynb