NaNs In The Output of SVD

Specific Domains Numerics
,
1

I am currently running numerical experiments which involve computing the SVD of a very large number of random matrices. Occasionally, Julia’s SVD function will return a Vt factor that contains a row of all NaNs. Here is some code which reproduces this behavior in Julia 1.8.5 and 1.8.1, using updated versions of LinearAlgebra and Random.

using LinearAlgebra
using Random

# constructing the matrix
rng = MersenneTwister(17, (0, 181935856, 181934854, 579))
D = Diagonal([ones(100); exp10.(range(0, -12, 50)); 1e-12*ones(100)])
Q = Matrix(qr(D*randn(rng, 250, 150)).Q)
A = Q'*D

svd(A).Vt # 139th row is all NaNs

This bug can be difficult to detect, since no error messages are thrown by the SVD to indicate an abnormal result. Running the SVD with alg = LinearAlgebra.QRIteration() fixes the issue, but I believe this algorithm is slower.

I hope this is useful information, and I am curious about the source of the error! Thank you for your time.

2 Likes
2

It may be a bug in the multithreading logic in OpenBLAS. I consistently get the NaN result with 6 or 8 threads, but no other number, so I suspect buffer mismanagement rather than a race condition.

I don’t see any NaNs when using MKL, which may also be faster.

2 Likes