Update 1 Per a request from Daniel Lichtblau at Wolfram Research, a minimal example has been uploaded to Google Drive as a pure-ASCII no-format "m"
-file named "oneMatrixThatFailsSVD_2017.m"
. The file is rather large (about 7.7 MBytes) because it encodes a convergence-failing 320x320
complex matrix as an integer byte-array of dimension 320x320x16
. This bit-perfect encoding "trick" is necessary because SVD convergence-failure seemingly is exquisitely sensitive to the least-significant bits of IEEE complex numbers (which are of course precisely the bits to which no well-conditioned SVD algorithm should be sensitive).
Update 2 A comment has been added to the end of the above-mentioned text file "oneMatrixThatFailsSVD_2017.m"
that provides the following (exact) pipe for flipping the least-significant-bit of complex matrices:
Flatten// (* complex square matrix -> complex list *)
ExportString[#,"Complex128"]&// (* complex list -> byte string *)
ExportString[#,"Base64"]&// (* byte string -> Base64 string *)
ByteArray[#]&//Normal// (* Base64 string -> integer list *)
Partition[#,16]&// (* integer list -> {...} *)
Partition[#,#//Length//Sqrt]&// (* {...} -> integer array *)
Map[(
{#[[1]]//If[#//EvenQ,#+1,#-1]&} ~ Join ~ (* flip LSB of real part *)
#[[2;;8]] ~ Join ~
{#[[9]]//If[#//EvenQ,#+1,#-1]&} ~ Join ~ (* flip LSB of imag part *)
#[[10;;16]]
)&,#,{2}]&//
Flatten// (* integer array -> integer list *)
ExportString[#,"Byte"]&// (* integer list -> byte string *)
ImportString[#,"Complex128"]&// (* byte string -> complex list *)
Partition[#,#//Length//Sqrt]&; (* complex list -> complex matrix *)
It turns out that flipping the least-significant bits of a convergence-failing input matrix does reproducibly eliminate the "SingularValueDecomposition::cflsvd:"
message (which is behavior that no well-conditioned SVD algorithm should exhibit).
Hopefully this fine-grained, exactly reproducible control of the "cflsvd"
SVD bug will make fixing it much easier in 2017 than back in 2005, when Mathematica tools like ByteArray[__]
and ExportString[_,"Complex128"]
were less-developed, such that the bug was more challenging to exhibit reproducibly and diagnose reliably.
The bug in a nutshell: for single-precision complex matrices, Mathematica's SingularValueDecomposition[] sporadically fails to converge.
The associated Mathematica-generated error message is:
SingularValueDecomposition::cflsvd: Machine-precision algorithm
failed to converge. Arbitrary-precision algorithm is called,
which is slower but more accurate.
This is a followup to a long-standing Mathematica bug report (specifically Wolfram Research bug report [TS 28968], submitted way back in 2005).
This tarball provides (in a folder named "SVDfailures_2017"
) 25 examples of matrices whose convergence fails under Mathematica 10.2.0 for Mac OS X x86. The same tarball provides (in a folder named "SVDfailures_2005
") matrices that fail under various versions of Mathematica dating back to 2005 (these files were provided with bug report [TS 28968]). The tarball is rather large (more than 100 MBytes) because it consists mostly of numerical matrices created with "DumpSave[__]"
).
A principal difference between the 2005 failures and the 2017 failures is that (at least some of) the matrices that failed outright back in 2005, now generate the (undocumented?) convergence-failure message "SingularValueDecomposition::cflsvd"
To anticipate some questions:
- The arbitrary-precision evaluation does yield a correct decomposition, at the expense of a runtime that is 500-1000X longer.
- The matrices that fail of convergence are (seemingly) unremarkable in respect to numerical condition and rank.
- There is no reason (known to me at least) why SVDs of
$\mathcal{O}(1)$-entry single-precision matrices should ever "fail to converge", and there is no linear algebra software other than Mathematica's (known to me) that exhibits a similar convergence failure.
My questions are:
- Does
"SingularValueDecomposition::cflsvd"
convergence-failure occur more generally, i.e., on systems other than Mathematica 10.2.0 for Mac OS X?
- What's the best way to report this bug (if it is a reproducible bug)?