[Numpy-discussion] Strange bug in SVD when numpy is installed with pip

[Darlan Cavalcante Moreira]

Some time ago I have reported a bug about the linalg.matrix_rank in
numpy for complex matrices. This was quickly fixed and to take the
advantage of the fix I'm now using numpy 1.9.0 installed through pip,
instead of the version from my system (Ubuntu 14.04, with numpy version
1.8.1).

However, I have now encountered a very strange bug in the SVD function,
but only when numpy is manually installed (with pip in my case). When I
calculate the SVD of complex matrices with more columns than rows the
last rows of the returned V_H matrix are all equal to zeros. This does
not happens for all shapes, but for the ones where this happens it will
always happen.

This can be reproduced with the code below. You can change the sizes of
M and N and it happens for other sizes where N > M, but now all of them.

--8<---------------cut here---------------start------------->8---
import numpy as np

M = 8                           # Number of rows
N = 12                          # Number of columns

# Calculate the SVD of a non-square complex random matrix
[U, S, V_H] = np.linalg.svd(np.random.randn(M, N) + 1j*np.random.randn(M, N), full_matrices=True)

# Calculate the norm of the submatrix formed by the last N-M rows
if np.linalg.norm(V_H[M-N:]) < 1e-30:
    print("Bug!")
else:
    print("No Bug")
# See the N-M rows. They are all equal to zeros
print(V_H[M-N:])
--8<---------------cut here---------------end--------------->8---

The original matrix can still be obtained from the decomposition, since
the zero rows correspond to zero singular values due to the fact that
the original matrix has more columns then rows. However, since the user
asked for 'full_matrices' here (the default) returning all zeroes for
these extra rows is not useful.

In order to isolate the bug I tried installing some different numpy
versions in different virtualenvs. I tried version 1.6, 1.7, 1.8.1 and
1.9 and the bug appears in all of then. Since it does not happen if I
use the version 1.8.1 installed through the Ubuntu package manager, I
imagine it is due to some issue when pip compiles numpy locally.

Note: If I run numpy.testing.test() all tests are OK for all numpy versions I
have tested. The only fail is a known fail and I get "OK (KNOWNFAIL=1)".


-- 
Darlan Cavalcante Moreira