book_kees
book_kees

Reputation: 347

Why is scipy's eigh returning unexpected negative eigenvalues?

When trying to calculate the eigenvalues of a particular matrix using scipy's eigh function, I am getting negative eigenvalues where I do not expect them. My code is as follows:

import numpy as np
from scipy import linalg

def get_symmetrized_eigens(array):
    a = np.dot(array.transpose(), array)  # Compute square matrix from the input

    evals, evecs = linalg.eigh(a)  # Compute eigenvalues, eigenvectors

    # Output to test for negative eigenvalues
    print('{0} eigenvalues computed'.format(evals.size))
    for x in evals:
        if x < 0:
            print(x)

    return evals, evecs

The code is supposed to take an arbitrarily-sized matrix (in ndarray format), generate a square matrix from it by left-multiplying its transpose, and calculate the eigenvalues and eigenvectors of the generated square matrix. Mathematically, the eigenvalues of a square matrix generated in this way must be greater than or equal to zero (given an m by n real-valued matrix A, let T be the tranpose of A. Then TA is a symmetric matrix that has real, non-negative eigenvalues). But when I input the following 35 by 40 matrix:

127 127 127 127 127 127 127 127 127 127 127 125 127 128 126 118 116 128 236 253 254 254 254 253 253 254 254 254 254 254 254 254 253 253 253 253 253 253 253 254
127 127 127 127 127 127 127 127 127 127 127 126 129 121 140 116 116 116 241 255 254 254 254 253 253 253 254 254 254 254 255 254 254 254 253 253 253 253 253 254
127 127 127 127 126 125 126 128 127 123 126 125 135 191 204 116 116 116 208 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255
127 127 127 126 123 127 129 125 124 129 126 169 233 252 173 116 116 116 173 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255
127 127 127 127 128 129 121 130 127 141 218 254 254 252 136 116 116 116 136 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255
127 127 127 124 125 124 131 125 174 243 252 253 252 237 116 116 116 116 116 239 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255
127 127 127 127 127 124 130 206 250 254 252 253 254 203 116 116 119 116 116 205 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255
127 127 127 126 126 138 227 253 254 254 252 254 253 169 116 116 173 116 116 170 255 255 255 255 255 255 255 255 255 255 255 255 249 245 255 255 255 255 255 255
122 117 117 119 123 152 165 174 183 190 199 208 217 132 116 117 237 117 116 132 217 209 200 192 183 175 165 156 146 137 127 118 158 238 255 255 255 255 255 255
127 127 118 116 116 116 116 116 116 116 116 116 116 116 116 147 255 147 116 116 116 116 116 116 116 116 116 116 116 116 129 209 255 255 255 255 255 255 255 255
125 126 148 170 116 116 116 116 116 116 116 116 116 116 116 185 255 184 116 116 116 116 116 116 116 116 116 116 116 173 245 255 255 255 255 255 255 255 255 255
125 143 239 252 221 139 116 116 116 116 153 189 180 172 163 226 255 226 163 172 180 189 153 116 116 116 116 139 222 255 255 255 255 255 255 255 255 255 255 255
137 238 255 255 255 251 187 119 116 116 116 164 242 255 255 255 255 255 255 255 242 164 116 116 116 119 187 251 255 255 255 255 255 255 255 255 255 255 255 255
227 254 252 255 255 255 255 232 150 116 116 116 125 201 254 255 255 255 254 201 125 116 116 116 150 232 255 255 255 255 255 255 255 255 255 255 255 255 255 255
254 252 252 255 255 255 255 255 254 199 124 116 116 116 220 255 255 255 220 116 116 116 124 199 254 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255
254 253 252 254 254 255 255 255 255 255 220 116 116 133 253 255 239 255 253 133 116 116 221 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255
254 252 253 254 255 255 255 255 255 255 181 116 116 189 255 190 117 190 255 189 116 116 181 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255
252 252 254 255 255 255 255 255 255 251 127 116 117 235 167 116 116 116 167 235 117 116 127 251 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255
253 254 253 255 255 255 255 255 255 209 116 116 143 146 116 116 116 116 116 146 142 116 116 209 252 252 252 254 253 254 253 251 255 252 251 255 255 255 255 255
254 252 254 255 255 255 255 255 255 153 116 116 116 116 116 120 177 120 116 116 116 116 116 153 253 252 254 253 253 254 252 253 251 251 253 255 255 255 255 255
254 251 253 255 255 255 255 255 235 116 116 116 116 116 129 225 255 225 129 116 116 116 116 116 235 254 251 252 254 253 254 253 254 251 253 255 255 255 255 255
253 253 254 255 255 255 255 255 184 116 116 116 116 144 240 255 255 255 240 143 116 116 116 116 181 253 253 254 254 254 253 252 253 253 250 255 255 255 255 255
250 254 253 255 255 255 255 252 129 116 116 116 163 249 255 255 255 255 255 247 162 116 116 116 121 176 241 251 254 253 254 254 253 253 254 255 255 255 255 255
253 254 252 255 255 255 255 211 116 116 116 187 255 255 255 255 255 255 255 253 252 186 116 116 116 126 151 229 250 252 252 254 254 251 254 255 255 255 255 255
254 253 255 255 255 255 255 156 116 121 209 255 255 255 255 255 255 255 255 253 253 255 207 121 116 123 125 124 190 244 253 252 253 254 253 255 255 255 255 255
254 254 255 255 255 255 238 117 132 228 255 255 255 255 255 255 255 255 255 251 253 253 254 225 121 116 128 128 123 163 225 253 254 254 253 255 255 255 255 255
255 255 255 255 255 255 186 147 242 255 255 255 255 255 255 255 255 255 255 254 254 254 254 254 149 121 124 125 129 125 133 194 252 253 253 254 255 253 254 253
255 255 255 255 255 252 181 250 255 255 255 255 255 255 255 255 255 255 255 254 254 254 254 254 152 127 124 127 126 130 127 125 167 228 253 250 253 253 253 253
255 255 255 255 255 252 255 255 255 255 255 255 255 255 255 255 255 255 255 254 254 254 254 254 152 126 126 130 123 124 125 131 124 137 199 253 251 254 252 254
255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 254 254 254 254 254 152 126 126 125 125 127 127 125 125 128 125 171 233 254 254 254
255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 254 254 254 254 254 152 126 126 127 127 126 125 124 125 126 127 125 142 209 254 253
255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 254 254 254 254 254 152 126 126 127 123 128 127 128 130 123 126 126 126 127 180 239
255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 254 254 254 254 254 152 126 126 124 129 123 126 128 124 127 127 127 123 128 127 146
255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 254 254 254 254 254 152 126 126 126 123 128 127 124 126 128 126 126 127 125 127 124
255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 254 254 254 254 254 152 126 126 127 127 127 127 127 127 127 127 126 127 125 126 128

I get the following output:

40 eigenvalues computed
-902.580414433
-829.138600111
-736.37232834
-649.343219906
-606.346570836
-542.284035259
-491.035503988
-433.836458775
-378.892895817
-349.983072146
-322.571278901
-272.422566695
-231.320418966
-215.794098609
-187.163923805
-171.74732025
-124.678786093
-83.6681245986
-36.464176635
-1.96859835522

So it seems to me that negative eigenvalues are being computed. I have searched far and wide, but similar questions (e.g. 1, 2) seem to deal with precision errors or different conditions on the matrix involved. Moreover, I had the same issue when using numpy's eigh function, as well as the more general eig function. Why is the calculation producing these negative eigenvalues, and how can I fix it?

Upvotes: 2

Views: 2465

Answers (1)

Śmigło
Śmigło

Reputation: 1037

Unfortunately this is because some lack of type control in numpy. REMEMBER: when something's wrong - check if you used something like ndarray.astype(np.int_) instead of ndarray.astype(np.int32).

You need to pass np.float64 or np.int32 matrix, not np.uint8 dtype here. Unfortunately such things matter in numpy and they can't be fixed and found easily... i've encountered it few times this year at work

so add a line:

def get_symmetrized_eigens(array):
    array = np.array(array, dtype=np.float64)
    a = np.dot(array.transpose(), array)  # Compute square matrix from the input

percusse wrote:

Having said that these routines are stress tested for decades now so it's very very unlikely that they choke on this array.

as you can see - it is likely :D

Upvotes: 1

Related Questions