Issue with Complex Matrices and Eigenvalues

[AjayKaladharan]

I am new to this package and have a problem similar to the one described in “05_matrix_operations.ipynb” (Example 5). My task involves a matrix with complex entries, where each entry has an associated error. I need to calculate the eigenvalues of this matrix, including the propagation of errors. The matrix is hermitian.

I attempted to define complex numbers with errors using “CObs,” but when I tried to compute the eigenvalues, I encountered the error: “‘CObs’ object has no attribute ‘cov_names’.” This seems to be a syntax issue, but I have been unable to resolve it.

I would greatly appreciate any guidance on how to correctly define a complex matrix with errors and perform matrix operations like eigenvalue decomposition.

[s-kuberski]

Hi,

thanks for your posting! Would it be possible for you to provide a minimal working example that leads to the error? This would make it a lot easier to debug.

[AjayKaladharan]

Hi,
Thank you for the quick response. Suppose I am interested in a 2×2 complex Hermitian matrix, where each entry has associated errors in both the real and imaginary parts. This is the code I used and the error I got. I would greatly appreciate any tips to resolve the issue.

%%%%%%%%%%%%%%%%%
import pyerrors as pe
import numpy as np

obs11 = pe.pseudo_Obs(4.1, 0.2, 'e1')
obs22 = pe.pseudo_Obs(1, 0.01, 'e2')
my_real_part = pe.pseudo_Obs(4.1, 0.2, 'e3')
my_imag_part = pe.pseudo_Obs(1, 0.01, 'e4')
my_imag_part0 = pe.pseudo_Obs(0, 0, 'e5')

my_cobs12 = pe.CObs(my_real_part, my_imag_part)
my_cobs21 = pe.CObs(my_real_part, -my_imag_part)
my_cobs11 = pe.CObs(obs11, my_imag_part0)
my_cobs22 = pe.CObs(obs11, my_imag_part0)
matrix = np.asarray([[my_cobs11,my_cobs12], [my_cobs12, my_cobs22]])
%%%%%%%%%%%%%%%%%%

The error message
10 e, v = pe.linalg.eigh(matrix)
11 pe.gm(e)
12 print('Eigenvalues:')

File /opt/anaconda3/lib/python3.11/site-packages/pyerrors/linalg.py:263, in eigh(obs, **kwargs)
261 def eigh(obs, **kwargs):
262 """Computes the eigenvalues and eigenvectors of a given hermitian matrix of Obs according to np.linalg.eigh."""
--> 263 w = derived_observable(lambda x, **kwargs: anp.linalg.eigh(x)[0], obs)
264 v = derived_observable(lambda x, **kwargs: anp.linalg.eigh(x)[1], obs)
265 return w, v

File /opt/anaconda3/lib/python3.11/site-packages/pyerrors/obs.py:1213, in derived_observable(func, data, array_mode, **kwargs)
1211 allcov = {}
1212 for o in raveled_data:
-> 1213 for name in o.cov_names:
1214 if name in allcov:
1215 if not np.allclose(allcov[name], o.covobs[name].cov):

AttributeError: 'CObs' object has no attribute 'cov_names'

[s-kuberski]

Hi,

I had a deep look at the code and tried a number of easy solutions, but so far it seems as if the problem is a bit more tricky than I expected it to be.

In case you are interested why this is the case: So far we circumvent the error propagation for complex values observables (CObs) in matrix operations by going to a real representation of the matrix. However, this does not work if one wants to solve the eigenproblem because the mapping is not just a linear algebra operation.
autograd seems to correctly propagate the derivatives with respect to the real and imaginary parts when using anp.linalg.eigh (as I have checked numerically). Based on this fact, one could implement the general error propagation for CObs, but this will be a complex addition to code and might take some time until one has derived the optimal approach.

Maybe I have overlooked something and @fjosw has apprach in mind.

If the functionality is needed, it could be worthwile to implement the feature. May I ask in what context you need to solve the eigenproblem? I never accountered this problem, I guess because correlation functions in quantum field theories with reflection positivity are either real or imaginary.

[fjosw]

I agree that this is a bit more tricky. I think one could maybe still go via a real representation and could make the error propagation for the real and imaginary parts of each eigenvalue separately before recombining them to CObs. That would produce complex conjugate pairs of eigenvalues and one would need to select the "correct" ones. I wrote a quick draft (without the selection):

import random
import autograd.numpy as anp
import pyerrors as pe

# Random complex matrix for testing
c_mat = np.array([pe.CObs(pe.pseudo_Obs(random.random(), random.random() * 0.1, "test"), pe.pseudo_Obs(random.random(), random.random() * 0.1, "test")) for _ in range(9)]).reshape(3, 3)

# Construct real block matrix representation
r_mat = np.vectorize(lambda x: x.real)(c_mat)
i_mat = np.vectorize(lambda x: x.imag)(c_mat)
block_matrix = np.block([[r_mat, -i_mat], [i_mat, r_mat]])


# Error propagation for real and imaginary parts of the eigenvalues
r_eig = pe.derived_observable(lambda x, **kwargs: anp.real(anp.linalg.eig(x)[0]), block_matrix)
i_eig = pe.derived_observable(lambda x, **kwargs: anp.imag(anp.linalg.eig(x)[0]), block_matrix)

# Combine eigenvalues to CObs
eigenvalues = []
for re, im in zip(r_eig, i_eig):
    eigenvalues.append(pe.CObs(re, im))
    eigenvalues[-1].gamma_method()
eigenvalues
>>> [CObs[(1.431(93)+1.663(53)j)],
 CObs[(1.431(93)-1.663(53)j)],
 CObs[(0.414(86)+0.630(57)j)],
 CObs[(0.414(86)-0.630(57)j)],
 CObs[(0.069(65)+0.105(61)j)],
 CObs[(0.069(65)-0.105(61)j)]]

# Verify eigenvalues
f_mat = np.vectorize(lambda x: float(x))(r_mat)+1j*np.vectorize(lambda x: float(x))(i_mat)
np.linalg.eig(f_mat).eigenvalues
>>> array([1.43118587+1.66349372j, 0.06898277-0.10543659j,
       0.41389627+0.63012354j])

Does that make sense?

[s-kuberski]

Oh wow, makes complete sense! I just noticed that there is a typo in the MWE that makes the matrix non-hermitian - I might have used that one when I tried (and failed with) the solution via the real representation, which is obviously inconsistent for eigh. I'll look into your solution, also for the eigenvectors and both eigh and eig.

[AjayKaladharan]

Hi, Thank you for your suggestion; it worked perfectly for the problem I am addressing. My work involves density matrices, where the matrix is complex Hermitian.

[s-kuberski]

Hi again,

this was a bit of a rabbit hole. In principle, it is possible to also implement a solution via the real representation that works for the eigenvectors. One has to cope with the possibility that different normalizations are used for the proper eigenproblem and the one with the real matrix, but this can be solved by projection.

However, in my tests with random (hermitian) matrices, I came across cases where the derivatives are all nan. A bit of googling shows that this is a known problem (also in jax), if the eigenvalues are (near-)degenerate. So it does not seem straight forward to implement a one-fits-all solution for the eigenvectors.

If we'd like to add the possibility to compute the eigenvalues of complex matrices (what was asked in this discussion), this could be implemented in pyerrors.linalg.eig (which is not a drop-in replacement for numpy.linalg.eig). Coping with the eigenvectors is more difficult and as long as no one needs it, we could put it on hold.