Real-valued results for sqrtm

[timholy]

The following matrix:

A = Float64[ 1421   52503   9933;
            52503 1942611 367521;
             9933  367521  69531]

poses problems for sqrtm. Currently, sqrtm gives a complex result, but that result is no more accurate than the real-valued version:

julia> M = sqrtm(A)
3x3 Array{Complex{Float64},2}:
 2.41272+3.26497e-30im  36.9632-9.27778e-19im  6.99303+4.90397e-18im
 36.9632-9.27778e-19im   1369.0+2.63639e-7im     259.0-1.39352e-6im 
 6.99303+4.90397e-18im    259.0-1.39352e-6im      49.0+7.36576e-6im 

julia> sumabs(A - M*M)
7.97695065557491e-9

julia> Mr = real(M)
3x3 Array{Float64,2}:
  2.41272    36.9632    6.99303
 36.9632   1369.0     259.0    
  6.99303   259.0      49.0    

julia> sumabs(A - Mr*Mr)
7.93329490988981e-9

The fundamental reason why it returns a complex-valued matrix seems to be this:

julia> F = eigfact(A)
Base.LinAlg.Eigen{Float64,Float64,Array{Float64,2},Array{Float64,1}}([-5.820766139459817e-11,1.9985905553312817,2.0135610014094408e6],3x3 Array{Float64,2}:
  6.54175e-13   0.999648   -0.0265466
 -0.185892     -0.0260839  -0.982224 
  0.98257      -0.0049348  -0.185826 )

julia> isposdef(F)
false

julia> F.values
3-element Array{Float64,1}:
 -5.82077e-11
  1.99859    
  2.01356e6  

On balance, I suspect we should make sqrtm more inclined to return a real-valued matrix. If others agree, I'm not sure where the right place to change this is, but my money is on adding a "sloppiness" here.

[timholy]

For completeness, I should have referenced relevant discussions in #4006.

[tkelman]

How to make the correct tolerancing decision on "sloppiness" in slightly negative eigenvalues is something that gives me pause. Would not want to do something like that without carefully thinking through the consequences.

edit: I'm wondering whether we might want to move towards user-provided positive definiteness annotations via the type system as a long-term direction for this.

[StefanKarpinski]

cc: @alanedelman

[timholy]

Agreed it's a slippery slope.

One data point: Matlab returns a real-valued matrix for this problem, but their eig also returns 2.5495e-10 for the smallest eigenvalue.

[andreasnoack]

If we consider your matrix a realization of random matrix, the smallest computed eigenvalues are symmetrically distributed around zero even though the true value is in fact zero. Therefore, the right choice of cutoff value or the sloppiness parameter would probably depend on the distribution of the entries and the true rank neither of which we know. It could be fun to take the statistical aspects of this seriously, but I guess a simplified approach is often preferred. E.g. it seems that MATLAB handles this in the same way as we are doing now. We use a different default eigensolvers, but we can reproduce MATLAB with

julia> LAPACK.syev!('V', 'U', copy(A))[1][1]
2.549678219012986e-10

so it is just chance that MATLAB happens to return a real result.

>> A = randi([1 100], 3,2);
>> sqrtm(A*A')

ans =

  42.4942 + 0.0000i  53.0067 - 0.0000i  40.1066 - 0.0000i
  53.0067 - 0.0000i  77.7233 + 0.0000i  35.7965 + 0.0000i
  40.1066 - 0.0000i  35.7965 + 0.0000i  55.3089 + 0.0000i

[timholy]

What about retrospectively deciding what to return? For example, does it make sense to take the real component and do something like sumabs(A - Mr*Mr), comparing that against a threshold? The problem I see with that test is that it's a second O(N^3) operation. It's one with a smaller coefficient (by about 4-fold in my tests), but it would still be a measurable burden.

[andreasnoack]

I think it is too costly to add such a test and effectively it must be very to setting small negative eigenvalues to zero which would be much cheaper. My feeling is that we should keep the present behavior unless we can come up with something more clever than the solutions mentioned so far. I you plan to roll your own positive semidefinite sqrtm then remember that you can ask for a specific part of the spectrum which will save you a little time, i.e. eigfact(Symmetric(A), 0.0, Inf).

[timholy]

Sounds reasonable to me, too.