# SVD: Better default to gesvd! instead of gesdd!?

**URL:** <https://discourse.julialang.org/t/svd-better-default-to-gesvd-instead-of-gesdd/20603>\
**Category:** Internals & Design\
**Tags:** linearalgebra\
**Created:** [February 9, 2019, 11:17am UTC](https://discourse.julialang.org/t/svd-better-default-to-gesvd-instead-of-gesdd/20603 "2019-02-09T11:17:53Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [February 9, 2019, 11:17am UTC](https://discourse.julialang.org/t/svd-better-default-to-gesvd-instead-of-gesdd/20603/1 "2019-02-09T11:17:53Z")

</div>

Currently, `svd` defaults to using `LinearAlgebra.LAPACK.gesdd!` under the hood. I [vaguely mentioned](https://discourse.julialang.org/t/call-lapacks-zheev-directly/7305/5) some time ago that `gesdd!` gives incorrect results in some cases which can be circumvented by explicitly calling `LinearAlgebra.LAPACK.gesvd!` instead. However, for large enough matrices, `gesvd!` is much slower than `gesdd!`, which [uses](https://www.netlib.org/lapack95/lug95/node324.html) a divide-and-conquer approach. Whether the inaccuracies (which can be huge) are due to a bug in the implementation or a shortcoming of this different method I don’t know.

In any case, it naturally raises the question which method we should default to. It seems to me that we should prefer accuracy over speed and default to `gesvd!` in `svd`. Of course, we could introduce a keyword flag `divide` or similar that allows one to switch to the faster method. Just getting wrong results without a warning or anything, as it is the case right now, seems bad.

Note that Octave, which was also defaulting to `gesdd` has recently made precisely this switch as discussed [in this issue](https://savannah.gnu.org/bugs/?55564).

Sort of taking the example from this thread, here is a MWE that demonstrates the issue:

```julia
using LinearAlgebra

function compan(c::AbstractVector{T}) where T
    N = length(c)
    A = diagm(-1 => ones(T, N-1))
    @inbounds for i in 1:N
        A[1,i] = - c[i]/c[1]
    end
    A
end

function svdvaldiff(N)
  A = compan(reverse( vcat(1, 1 ./ cumprod(1:N, dims=1)) ));
  svals_gesdd = LinearAlgebra.LAPACK.gesdd!('A', copy(A))[2];
  svals_gesvd = LinearAlgebra.LAPACK.gesvd!('A', 'A', copy(A))[2];

  return maximum(svals_gesdd - svals_gesvd)
end

function gesdd_gesvd_comparison(Nmax=40)
    for N in 1:Nmax
        println("N = $(N), max diff svdvals: ", svdvaldiff(N)) 
    end
end

```

Which gives the following output:

```julia
julia> gesdd_gesvd_comparison()
N = 1, max diff svdvals: 0.0
N = 2, max diff svdvals: 0.0
N = 3, max diff svdvals: 0.0
N = 4, max diff svdvals: 0.0
N = 5, max diff svdvals: 0.0
N = 6, max diff svdvals: 0.0
N = 7, max diff svdvals: 0.0
N = 8, max diff svdvals: 0.0
N = 9, max diff svdvals: 0.0
N = 10, max diff svdvals: 0.0
N = 11, max diff svdvals: 0.0
N = 12, max diff svdvals: 0.0
N = 13, max diff svdvals: 0.0
N = 14, max diff svdvals: 0.0
N = 15, max diff svdvals: 0.0
N = 16, max diff svdvals: 0.0
N = 17, max diff svdvals: 0.0
N = 18, max diff svdvals: 0.0
N = 19, max diff svdvals: 0.0
N = 20, max diff svdvals: 0.0
N = 21, max diff svdvals: 0.0
N = 22, max diff svdvals: 0.0
N = 23, max diff svdvals: 0.0
N = 24, max diff svdvals: 0.0
N = 25, max diff svdvals: 749.7516745517077
N = 26, max diff svdvals: 166.767736950586
N = 27, max diff svdvals: 584.3097624800283
N = 28, max diff svdvals: 635.9955873150014
N = 29, max diff svdvals: 752.0372910883041
N = 30, max diff svdvals: 934.3056077045787
N = 31, max diff svdvals: 532.6321090336445
N = 32, max diff svdvals: 644.2861698543643
N = 33, max diff svdvals: 362.0563181932318
N = 34, max diff svdvals: 524.8519471908377
N = 35, max diff svdvals: 681.9577012165146
N = 36, max diff svdvals: 959.8179240905338
N = 37, max diff svdvals: 116.35110997021529
N = 38, max diff svdvals: 510.2091542800206
N = 39, max diff svdvals: 245.12683184534578
N = 40, max diff svdvals: 6.869915826281117

```

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [February 9, 2019, 11:38am UTC](https://discourse.julialang.org/t/svd-better-default-to-gesvd-instead-of-gesdd/20603/2 "2019-02-09T11:38:09Z")

</div>

> [@carstenbauer](#):
>
> we should prefer accuracy over speed and default to `gesvd!` in `svd` .

I agree. The Octave discussion you linked mentions a LAPACK issue,

> <https://github.com/Reference-LAPACK/lapack/issues/316>
>
> Last night I found what appears to be a bad bug when computing a full SVD using …xGESDD. On the particular matrix, companion\_demo(n) from EigTool with n\>=26, the computed singular values are totally different than those when also requesting the singular vectors. When only the singular values are requested, both xGESDD and xGESVD seem to work fine, regardless of n. But when n\>=26 and singular vectors are requested, xGESDD fails.
> 
> I actually ran into this problem using MATLAB and wrote a demo to replicate the issue. Thanks to Martin Köhler for quickly translating my demo to Fortran. Both demos, with documentation, are provided. To run the MATLAB demo, you will also need to install EigTool to generate the problematic companion matrix. Just drop the .txt extensions on the Fortran and MATLAB demos (added so I could upload them here).
> 
> Many thanks,
> Tim
> 
> \[README.txt\](https://github.com/Reference-LAPACK/lapack/files/2788012/README.txt)
> 
> \[svdbug.f90.txt\](https://github.com/Reference-LAPACK/lapack/files/2788020/svdbug.f90.txt)
> 
> \[svdbug.m.txt\](https://github.com/Reference-LAPACK/lapack/files/2788021/svdbug.m.txt)

so they are working on the problem, and may solve/clarify it by the time of the next release **EDIT** _of Julia_. Perhaps you could open an issue (for Julia) and xref the above one.

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [February 9, 2019, 1:46pm UTC](https://discourse.julialang.org/t/svd-better-default-to-gesvd-instead-of-gesdd/20603/3 "2019-02-09T13:46:25Z")

</div>

I should probably clarify one more thing. The bug (or whatever) in `gesdd!` seems to only be present in full mode (`'A'`), as indicated in the octave/LAPACK issues mentioned above. If one replaces the `gesdd!` call in the MWE by

```julia
svals_gesdd = LinearAlgebra.LAPACK.gesdd!('N', copy(A))[2]; # 'N' instead of 'A'

```

one gets

```julia
N = 1, max diff svdvals: 0.0
N = 2, max diff svdvals: 1.1102230246251565e-16
N = 3, max diff svdvals: 0.0
N = 4, max diff svdvals: 7.105427357601002e-15
N = 5, max diff svdvals: 2.842170943040401e-14
N = 6, max diff svdvals: 2.220446049250313e-16
N = 7, max diff svdvals: 0.0
N = 8, max diff svdvals: 2.220446049250313e-16
N = 9, max diff svdvals: 1.1102230246251565e-16
N = 10, max diff svdvals: 9.313225746154785e-10
N = 11, max diff svdvals: 0.0
N = 12, max diff svdvals: 0.0
N = 13, max diff svdvals: 5.551115123125783e-16
N = 14, max diff svdvals: 2.220446049250313e-16
N = 15, max diff svdvals: 0.0
N = 16, max diff svdvals: 0.0
N = 17, max diff svdvals: 2.220446049250313e-16
N = 18, max diff svdvals: 2.220446049250313e-16
N = 19, max diff svdvals: 32.0
N = 20, max diff svdvals: 1.1368683772161603e-13
N = 21, max diff svdvals: 1.1368683772161603e-13
N = 22, max diff svdvals: 2.220446049250313e-16
N = 23, max diff svdvals: 2.220446049250313e-16
N = 24, max diff svdvals: 2.220446049250313e-16
N = 25, max diff svdvals: 5.684341886080802e-14
N = 26, max diff svdvals: 4.336808689942018e-19
N = 27, max diff svdvals: 5.684341886080802e-14
N = 28, max diff svdvals: 1.1102230246251565e-16
N = 29, max diff svdvals: 2.2737367544323206e-13
N = 30, max diff svdvals: 1.1368683772161603e-13
N = 31, max diff svdvals: 5.684341886080802e-14
N = 32, max diff svdvals: 2.220446049250313e-16
N = 33, max diff svdvals: 1.1368683772161603e-13
N = 34, max diff svdvals: 1.1102230246251565e-16
N = 35, max diff svdvals: 9.094947017729282e-13
N = 36, max diff svdvals: 1.1102230246251565e-16
N = 37, max diff svdvals: 2.842170943040401e-14
N = 38, max diff svdvals: 1.1102230246251565e-16
N = 39, max diff svdvals: 1.4210854715202004e-14
N = 40, max diff svdvals: 2.220446049250313e-16

```

Hence, only Julia’s `svd(A, full=true)` is effected.

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [February 9, 2019, 2:02pm UTC](https://discourse.julialang.org/t/svd-better-default-to-gesvd-instead-of-gesdd/20603/4 "2019-02-09T14:02:03Z")

</div>

Done, see [`svd` gives wrong results (upstream LAPACK bug) · Issue #31023 · JuliaLang/julia · GitHub](https://github.com/JuliaLang/julia/issues/31023).

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [February 9, 2019, 2:59pm UTC](https://discourse.julialang.org/t/svd-better-default-to-gesvd-instead-of-gesdd/20603/5 "2019-02-09T14:59:29Z")

</div>

That matrix appears to be pathologically badly conditioned: it has some entries around 1e16, and you’re looking at singular values of order unity. This means that a _relative_ perturbation of the matrix of size 1e-16, as might result from machine precision, will appreciably change singular values around unity. So I’m not sure if this should be called a bug rather than the normal effects of machine precision; sacrificing performance in the non-pathological case seems premature.

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [February 9, 2019, 3:03pm UTC](https://discourse.julialang.org/t/svd-better-default-to-gesvd-instead-of-gesdd/20603/6 "2019-02-09T15:03:02Z")

</div>

Also I can’t seem to reproduce in pure julia, did I do something wrong?

```julia
N = 37; norm(svd(compan(reverse( vcat(1, 1 ./ cumprod(1:N, dims=1)) )),full=true).S - svd(compan(reverse( vcat(1, 1 ./ cumprod(1:N, dims=1)) ))).S)

```

gets me 0

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [February 9, 2019, 3:12pm UTC](https://discourse.julialang.org/t/svd-better-default-to-gesvd-instead-of-gesdd/20603/7 "2019-02-09T15:12:51Z")

</div>

I’m on mobile, but I think you shouldn’t calculate the norm but compare singular values individually.

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [February 9, 2019, 3:16pm UTC](https://discourse.julialang.org/t/svd-better-default-to-gesvd-instead-of-gesdd/20603/8 "2019-02-09T15:16:16Z")

</div>

I took the matrix from the octave thread. However, I’ve faced this issue in a concrete physics simulation before. Badly conditioned matrices do show up in real world applications. In fact, I use svd decompositions in these scenarios to avoid mixing of scales in matrix products and in matrix inversion.

---

<div class="post-metadata">

**Author:** ![Ralph\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ralph_smith/32/10344_2.png) [@Ralph\_Smith](https://discourse.julialang.org/u/Ralph_Smith)\
**Post date:** [February 9, 2019, 3:25pm UTC](https://discourse.julialang.org/t/svd-better-default-to-gesvd-instead-of-gesdd/20603/9 "2019-02-09T15:25:56Z")

</div>

These matrices are indeed nasty. **Neither** LAPACK routine gives an accurate decomposition for N\>25 because of rounding errors for `Float64` . Fortunately in Julia we have [GenericSVD](https://github.com/JuliaMath/GenericSVD.jl) which does a good job with `BigFloat` matrices.

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [February 9, 2019, 3:27pm UTC](https://discourse.julialang.org/t/svd-better-default-to-gesvd-instead-of-gesdd/20603/10 "2019-02-09T15:27:49Z")

</div>

They are identical. Norm(a-b) =0 means a==b, no trickiness here.

Yes, badly conditioned matrices are a fact of life, and one has to deal with them to avoid loss of precision. But there’s not much more lapack can do here I think - different methods have slightly different results, both compatible with the error bars expected. You’re likely to see similar differences on different architectures, threading, etc.

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [February 9, 2019, 3:57pm UTC](https://discourse.julialang.org/t/svd-better-default-to-gesvd-instead-of-gesdd/20603/11 "2019-02-09T15:57:22Z")

</div>

Of course, my bad…

Turns out Julia is inconsistent: `norm(svdvals(A) - svd(A).S) != 0`, independent of the value of `full`.

However, if I look at [the implementation](https://github.com/JuliaLang/julia/blob/master/stdlib/LinearAlgebra/src/svd.jl#L63) full should switch between `'A'` and `'N'`. Don’t know what’s going on here.

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [February 9, 2019, 4:08pm UTC](https://discourse.julialang.org/t/svd-better-default-to-gesvd-instead-of-gesdd/20603/12 "2019-02-09T16:08:58Z")

</div>

I’m far from an expert on this, so someone more qualified than me should judge whether this is a bug or not. In my physics simulation switching from `gesdd!` to `gesvd!` made all the difference. Before this change everything was wrong, afterwards the results were correct. So it’s hard for me to be satisfied with “neither routine works”/“they are equally bad”.

In this direction, how would you explain/comment the following observation:

```julia
N = 37
A = compan(reverse( vcat(1, 1 ./ cumprod(1:N, dims=1)) ))

using GenericSVD
B = BigFloat.(A)
vals = svd(B).S

display("gesdd:")
display(Float64.(LinearAlgebra.LAPACK.gesdd!('A', A)[2] - vals))
println()

display("gesvd: ")
display(Float64.(LinearAlgebra.LAPACK.gesvd!('A', 'A', A)[2] - vals))
println()

```

which gives

```julia
"gesdd:"
38-element Array{Float64,1}:
 -545.9257605455043 
  180.8251664583235 
  133.97039557171414
  133.97039557171414
  133.97039557171414
  133.97039557171414
  133.97039557171414
  133.97039557171414
  133.97039557171414
  133.97039557171414
  133.97039557171414
  133.97039557171414
  133.97039557171414
    ⋮               
  133.97039557171414
  133.97039557171414
  133.97039557171414
  133.97039557171414
  133.97039557171414
  133.97039557171414
  133.97039557171414
  133.97039557171414
  133.97039557171414
  133.97039557171414
  133.97039557171414
   81.05435621882137
"gesvd: "
38-element Array{Float64,1}:
 -545.9257605455043    
   95.97867447307135   
    5.3787786348333455 
    5.156126423112092  
    5.006756486183942  
    4.871019085393854  
    4.461591733737439  
    4.181477907103048  
    3.980235090671698  
    3.84337889798058   
    3.7403253010761617 
    3.5582886967320233 
    3.3537939211414782 
    ⋮                  
    0.8753460484851945 
    0.7662725689736753 
    0.664655536722611  
    0.5026707900598488 
    0.40756319405669306
    0.3737408827688582 
   -0.0506366143480812 
   -0.13352353403151718
   -0.2954480014961395 
   -0.5336271002375245 
   -0.77719691671636   
   -0.655727106867289  

```

The latter looks much better to me.

---

<div class="post-metadata">

**Author:** ![Ralph\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ralph_smith/32/10344_2.png) [@Ralph\_Smith](https://discourse.julialang.org/u/Ralph_Smith)\
**Post date:** [February 9, 2019, 5:00pm UTC](https://discourse.julialang.org/t/svd-better-default-to-gesvd-instead-of-gesdd/20603/13 "2019-02-09T17:00:00Z")

</div>

The [error bounds for SVD](http://www.netlib.org/lapack/lug/node96.html) imply that the singular values are only good to within about `±350` for N=37 in `Float64`, so both routines are actually performing as advertised in that regard. (One might quibble about the largest s.v. in both cases on your machine, but it’s still close.) In this particular case, one is lucky that many singular values are closer for `gesvd!` (in fact on my machine they are even closer than your results, consistent with @antoine-levitt’s point above). But note that the relative errors for the second-largest and smallest s.v. are pretty bad for both, so common applications of the SVD would be misleading here. This just isn’t a good problem for `Float64`.

I don’t mean to assert that `gesdd!` and `gesvd!` are generally equally reliable, just that this example is not conclusive. It sounds as if you have some interesting matrices in your applications; if you could post some examples they might help to elucidate the difference.

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [February 9, 2019, 5:09pm UTC](https://discourse.julialang.org/t/svd-better-default-to-gesvd-instead-of-gesdd/20603/14 "2019-02-09T17:09:36Z")

</div>

In fact a similar case is discussed explicitly in [http://www.netlib.org/lapack/lug/node97.html:](http://www.netlib.org/lapack/lug/node97.html:) “Thus large singular values (those near sigma\_1) are computed to high relative accuracy and small ones may not be.”. The surprise here for me is actually that gesvd is performing better than expected, which is good but probably shouldn’t be relied upon.

What’s your exact application, ie what do you need to do to the matrix?

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [February 10, 2019, 9:54am UTC](https://discourse.julialang.org/t/svd-better-default-to-gesvd-instead-of-gesdd/20603/15 "2019-02-10T09:54:10Z")

</div>

> [@carstenbauer](#):
>
> we should prefer accuracy over speed

This is inconsistent with `\ `, which uses LU instead of QR because it prefers speed over accuracy. So probably the default should remain the faster one.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [February 10, 2019, 10:06am UTC](https://discourse.julialang.org/t/svd-better-default-to-gesvd-instead-of-gesdd/20603/16 "2019-02-10T10:06:46Z")

</div>

I think that in general these decisions should be more nuanced, and depend on the degree of inaccuracy and the speed gain.

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [February 10, 2019, 4:52pm UTC](https://discourse.julialang.org/t/svd-better-default-to-gesvd-instead-of-gesdd/20603/17 "2019-02-10T16:52:29Z")

</div>

Ok, I managed to revive my 2 year old code and extract a MWE that displays a qualitative difference between `gesdd!` and `gesvd!` when used to stabilize a matrix multiplication.

Given a matrix `B`, the task is to calculate the matrix product chain `B * B * B * B * ...`. (Eventually I want to invert this chain and calculate the determinant of the inverse, but this doesn’t matter here I guess).

If we perform the multiplication naively (left plot) we’ll loose the lowest singular values due to mixing of scales. Note that the finite precision artifacts (at around 60 multiplications) kick in when the lowest singular values reach roughly the floating point precision limitation.

To avoid this issue, we stabilize the procedure by SVD decomposing in every iteration, separating the scales. However, only when we use `gesvd!` (right plot) do we see the expected result. When using `gesdd!` (middle) instead the problem still persists.

 ![product_chain_stable](https://global.discourse-cdn.com/julialang/original/3X/c/9/c94f8a19813bd3d3d1df283dd433e890d1ac614d.png)

MWE:

```julia
using PyCall, PyPlot, LinearAlgebra
@pyimport matplotlib.ticker as ticker

B = [0.9631774046507462 0.048118777932499635 0.002403935950466714 0.04811877793249931 0.09599796038960963 0.004795902100335193 0.00023959547538980293 0.004795902100335332 0.00956792420011803 0.00047799794475834705 2.3880000553206727e-5 0.00047799794475778125 0.0959979603896096 0.004795902100335486 0.0002395954753902608 0.004795902100335878 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.048118777932499635 0.9631774046507456 0.04811877793249987 0.002403935950467196 0.004795902100335861 0.09599796038960962 0.0047959021003353515 0.0002395954753906275 0.0004779979447579591 0.00956792420011812 0.0004779979447577643 2.388000055312064e-5 0.0047959021003361885 0.09599796038960963 0.0047959021003355865 0.00023959547538973634 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.002403935950466714 0.04811877793249987 0.9631774046507452 0.048118777932499295 0.00023959547538959772 0.0047959021003353254 0.0959979603896097 0.004795902100335392 2.388000055344459e-5 0.0004779979447581358 0.009567924200118488 0.0004779979447580419 0.0002395954753897907 0.004795902100335456 0.09599796038961016 0.004795902100335967 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.04811877793249931 0.002403935950467196 0.048118777932499295 0.9631774046507461 0.0047959021003360055 0.0002395954753902414 0.004795902100335805 0.095997960389609 0.0004779979447578188 2.3880000553252446e-5 0.00047799794475820914 0.00956792420011871 0.0047959021003361755 0.00023959547539000118 0.0047959021003356125 0.09599796038960898 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.09599796038960963 0.004795902100335861 0.00023959547538959772 0.0047959021003360055 0.9631774046507459 0.04811877793249995 0.0024039359504669402 0.04811877793249928 0.09599796038960942 0.004795902100335522 0.00023959547539085065 0.004795902100336073 0.009567924200118628 0.0004779979447581598 2.388000055308229e-5 0.0004779979447580583 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.004795902100335193 0.09599796038960962 0.0047959021003353254 0.0002395954753902414 0.04811877793249995 0.9631774046507461 0.04811877793250019 0.0024039359504671003 0.0047959021003356785 0.09599796038960914 0.004795902100335187 0.00023959547539012586 0.000477997944758397 0.009567924200118712 0.0004779979447579929 2.388000055378612e-5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.00023959547538980293 0.0047959021003353515 0.0959979603896097 0.004795902100335805 0.0024039359504669402 0.04811877793250019 0.9631774046507459 0.048118777932499454 0.00023959547538995215 0.004795902100335738 0.09599796038960945 0.004795902100335804 2.388000055351635e-5 0.00047799794475794286 0.009567924200118543 0.00047799794475805626 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.004795902100335332 0.0002395954753906275 0.004795902100335392 0.095997960389609 0.04811877793249928 0.0024039359504671003 0.048118777932499454 0.9631774046507461 0.00479590210033561 0.0002395954753902095 0.0047959021003353905 0.0959979603896089 0.0004779979447580881 2.3880000552118516e-5 0.0004779979447578195 0.009567924200117978 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.00956792420011803 0.0004779979447579591 2.388000055344459e-5 0.0004779979447578188 0.09599796038960942 0.0047959021003356785 0.00023959547538995215 0.00479590210033561 0.9631774046507453 0.048118777932500585 0.002403935950466673 0.0481187779325001 0.09599796038960902 0.004795902100335514 0.00023959547538961477 0.004795902100335855 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.00047799794475834705 0.00956792420011812 0.0004779979447581358 2.3880000553252446e-5 0.004795902100335522 0.09599796038960914 0.004795902100335738 0.0002395954753902095 0.048118777932500585 0.9631774046507456 0.04811877793249896 0.002403935950467324 0.004795902100335815 0.09599796038960987 0.004795902100335308 0.00023959547538959935 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 2.3880000553206727e-5 0.0004779979447577643 0.009567924200118488 0.00047799794475820914 0.00023959547539085065 0.004795902100335187 0.09599796038960945 0.0047959021003353905 0.002403935950466673 0.04811877793249896 0.9631774046507457 0.048118777932499086 0.00023959547539017492 0.0047959021003355804 0.09599796038960905 0.004795902100335983 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.00047799794475778125 2.388000055312064e-5 0.0004779979447580419 0.00956792420011871 0.004795902100336073 0.00023959547539012586 0.004795902100335804 0.0959979603896089 0.0481187779325001 0.002403935950467324 0.048118777932499086 0.963177404650746 0.004795902100335326 0.00023959547538994093 0.004795902100335743 0.09599796038960902 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0959979603896096 0.0047959021003361885 0.0002395954753897907 0.0047959021003361755 0.009567924200118628 0.000477997944758397 2.388000055351635e-5 0.0004779979447580881 0.09599796038960902 0.004795902100335815 0.00023959547539017492 0.004795902100335326 0.9631774046507464 0.048118777932499905 0.0024039359504668153 0.04811877793249932 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.004795902100335486 0.09599796038960963 0.004795902100335456 0.00023959547539000118 0.0004779979447581598 0.009567924200118712 0.00047799794475794286 2.3880000552118516e-5 0.004795902100335514 0.09599796038960987 0.0047959021003355804 0.00023959547538994093 0.048118777932499905 0.963177404650745 0.048118777932499884 0.002403935950466461 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0002395954753902608 0.0047959021003355865 0.09599796038961016 0.0047959021003356125 2.388000055308229e-5 0.0004779979447579929 0.009567924200118543 0.0004779979447578195 0.00023959547538961477 0.004795902100335308 0.09599796038960905 0.004795902100335743 0.0024039359504668153 0.048118777932499884 0.9631774046507459 0.048118777932499315 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.004795902100335878 0.00023959547538973634 0.004795902100335967 0.09599796038960898 0.0004779979447580583 2.388000055378612e-5 0.00047799794475805626 0.009567924200117978 0.004795902100335855 0.00023959547538959935 0.004795902100335983 0.09599796038960902 0.04811877793249932 0.002403935950466461 0.048118777932499315 0.9631774046507473 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.9631774046507459 -0.0959979603896101 0.009567924200118063 -0.09599796038960919 -0.04811877793249977 0.0047959021003357106 -0.0004779979447580722 0.004795902100335954 0.0024039359504662056 -0.00023959547538962708 2.38800005536094e-5 -0.00023959547539008502 -0.04811877793249954 0.004795902100335876 -0.0004779979447581868 0.004795902100335842; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0959979603896101 0.9631774046507463 -0.09599796038960923 0.00956792420011795 0.00479590210033555 -0.04811877793250019 0.004795902100335323 -0.0004779979447578979 -0.00023959547539012795 0.0024039359504670526 -0.00023959547538991325 2.3880000553715824e-5 0.004795902100335467 -0.04811877793249935 0.004795902100335896 -0.00047799794475811004; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.009567924200118063 -0.09599796038960923 0.963177404650746 -0.09599796038960962 -0.0004779979447578395 0.004795902100335537 -0.04811877793249961 0.004795902100335916 2.3880000553300117e-5 -0.00023959547539028687 0.0024039359504665187 -0.00023959547538977485 -0.00047799794475783845 0.004795902100335813 -0.04811877793249962 0.004795902100335811; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.09599796038960919 0.00956792420011795 -0.09599796038960962 0.9631774046507457 0.00479590210033584 -0.00047799794475798536 0.004795902100336278 -0.048118777932499524 -0.00023959547538984473 2.388000055350761e-5 -0.0002395954753899119 0.0024039359504664567 0.004795902100335905 -0.00047799794475835626 0.004795902100336193 -0.048118777932499614; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.04811877793249977 0.00479590210033555 -0.0004779979447578395 0.00479590210033584 0.9631774046507465 -0.09599796038960896 0.009567924200117685 -0.09599796038960907 -0.048118777932499954 0.0047959021003352465 -0.00047799794475821906 0.004795902100336257 0.0024039359504666635 -0.00023959547539018227 2.388000055351761e-5 -0.0002395954753895482; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0047959021003357106 -0.04811877793250019 0.004795902100335537 -0.00047799794475798536 -0.09599796038960896 0.9631774046507452 -0.09599796038960884 0.00956792420011878 0.004795902100335861 -0.04811877793250033 0.004795902100335711 -0.00047799794475756154 -0.00023959547538969806 0.002403935950466428 -0.00023959547538983226 2.3880000553934592e-5; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0004779979447580722 0.004795902100335323 -0.04811877793249961 0.004795902100336278 0.009567924200117685 -0.09599796038960884 0.9631774046507464 -0.09599796038960885 -0.0004779979447580733 0.004795902100335925 -0.04811877793250011 0.004795902100335748 2.388000055362013e-5 -0.00023959547539026475 0.0024039359504664805 -0.00023959547538962662; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.004795902100335954 -0.0004779979447578979 0.004795902100335916 -0.048118777932499524 -0.09599796038960907 0.00956792420011878 -0.09599796038960885 0.963177404650745 0.004795902100335986 -0.0004779979447581103 0.004795902100335992 -0.04811877793249953 -0.00023959547539023258 2.3880000553186046e-5 -0.0002395954753903041 0.002403935950466679; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0024039359504662056 -0.00023959547539012795 2.3880000553300117e-5 -0.00023959547538984473 -0.048118777932499954 0.004795902100335861 -0.0004779979447580733 0.004795902100335986 0.9631774046507452 -0.09599796038960969 0.009567924200118537 -0.09599796038960912 -0.04811877793249956 0.00479590210033607 -0.0004779979447581861 0.00479590210033581; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.00023959547538962708 0.0024039359504670526 -0.00023959547539028687 2.388000055350761e-5 0.0047959021003352465 -0.04811877793250033 0.004795902100335925 -0.0004779979447581103 -0.09599796038960969 0.9631774046507461 -0.0959979603896091 0.009567924200117853 0.004795902100335603 -0.04811877793249963 0.00479590210033595 -0.00047799794475808386; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 2.38800005536094e-5 -0.00023959547538991325 0.0024039359504665187 -0.0002395954753899119 -0.00047799794475821906 0.004795902100335711 -0.04811877793250011 0.004795902100335992 0.009567924200118537 -0.0959979603896091 0.9631774046507453 -0.09599796038960987 -0.0004779979447580813 0.004795902100335885 -0.048118777932499204 0.004795902100335838; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.00023959547539008502 2.3880000553715824e-5 -0.00023959547538977485 0.0024039359504664567 0.004795902100336257 -0.00047799794475756154 0.004795902100335748 -0.04811877793249953 -0.09599796038960912 0.009567924200117853 -0.09599796038960987 0.9631774046507457 0.004795902100335232 -0.00047799794475815053 0.004795902100336965 -0.04811877793249964; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.04811877793249954 0.004795902100335467 -0.00047799794475783845 0.004795902100335905 0.0024039359504666635 -0.00023959547538969806 2.388000055362013e-5 -0.00023959547539023258 -0.04811877793249956 0.004795902100335603 -0.0004779979447580813 0.004795902100335232 0.9631774046507473 -0.09599796038960959 0.009567924200117964 -0.0959979603896097; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.004795902100335876 -0.04811877793249935 0.004795902100335813 -0.00047799794475835626 -0.00023959547539018227 0.002403935950466428 -0.00023959547539026475 2.3880000553186046e-5 0.00479590210033607 -0.04811877793249963 0.004795902100335885 -0.00047799794475815053 -0.09599796038960959 0.9631774046507461 -0.09599796038960948 0.00956792420011772; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0004779979447581868 0.004795902100335896 -0.04811877793249962 0.004795902100336193 2.388000055351761e-5 -0.00023959547538983226 0.0024039359504664805 -0.0002395954753903041 -0.0004779979447581861 0.00479590210033595 -0.048118777932499204 0.004795902100336965 0.009567924200117964 -0.09599796038960948 0.9631774046507465 -0.09599796038960973; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.004795902100335842 -0.00047799794475811004 0.004795902100335811 -0.048118777932499614 -0.0002395954753895482 2.3880000553934592e-5 -0.00023959547538962662 0.002403935950466679 0.00479590210033581 -0.00047799794475808386 0.004795902100335838 -0.04811877793249964 -0.0959979603896097 0.00956792420011772 -0.09599796038960973 0.9631774046507472]

function calc_product_chain(B, N)
    bs = size(B,1)
    R = Matrix{Float64}(I, bs, bs)
    U = Matrix{Float64}(I, bs, bs)
    S = ones(bs)
    Vt = Matrix{Float64}(I, bs, bs)

    svs = zeros(bs,N)
    svc = 1
    for k in 1:N
        R = B * R
        U, S, Vt = LinearAlgebra.LAPACK.gesvd!('A', 'A', copy(R))
        svs[:,svc] = log.(S)
        svc += 1
    end
    return (R, svs)
end

function calc_product_chain_stabilized(B, N, svdmethod)
    bs = size(B,1)

    U = Matrix{Float64}(I, bs, bs)
    S = ones(bs)
    Vt = Matrix{Float64}(I, bs, bs)
    Vtnew = Matrix{Float64}(I, bs, bs)

    svs = zeros(bs,N)
    svc = 1
    for k in 1:N
        U = B * U
        U *= Diagonal(S)
        U, S, Vtnew = svdmethod(U)
        Vt = Vtnew * Vt
        svs[:,svc] = log.(S)
        svc += 1
    end
    return (U,S,Vt,svs)
end

function plot_product_chain_svs(B, N)
  T = calc_product_chain(B, N)
  svs = T[2]
  fig = figure(figsize=(20,7))
  b = fig[:add_subplot](131)
  b[:yaxis][:set_major_locator](ticker.MaxNLocator(symmetric=true))
  b[:plot](copy(svs[:,:]'))
  b[:set_xlabel]("# multiplications")
  b[:set_ylabel]("log singular values of matrix product chain")
  b[:set_title]("naive")
    
  T = calc_product_chain_stabilized(B, N, x -> LinearAlgebra.LAPACK.gesdd!('A', x))
  svs = T[4]
  c = fig[:add_subplot](132, sharey=b, sharex=b)
  c[:plot](copy(svs[:,:]'))
  setp(c[:get_yticklabels](), visible=false)
  c[:set_xlabel]("# multiplications")
  c[:set_title]("stabilized, gesdd!")
    
  T = calc_product_chain_stabilized(B, N, x -> LinearAlgebra.LAPACK.gesvd!('A', 'A', x))
  svs = T[4]
  c = fig[:add_subplot](133, sharey=b, sharex=b)
  c[:plot](copy(svs[:,:]'))
  setp(c[:get_yticklabels](), visible=false)
  c[:set_xlabel]("# multiplications")
  c[:set_title]("stabilized, gesvd!")

  # tight_layout()
  subplots_adjust(wspace = 0.0)
  show()
  println(maximum(svs))
  nothing
end

plot_product_chain_svs(B, 150)

```

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [February 10, 2019, 5:42pm UTC](https://discourse.julialang.org/t/svd-better-default-to-gesvd-instead-of-gesdd/20603/18 "2019-02-10T17:42:51Z")

</div>

Your matrix is symmetric and you’re interested in repeated multiplications, so don’t you want an eigenvalue decomposition rather than an SVD? I’m not sure I understand your code but it looks like you’re assuming that U==Vt, which is not the case e.g. for negative eigenvalues.

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [February 10, 2019, 6:07pm UTC](https://discourse.julialang.org/t/svd-better-default-to-gesvd-instead-of-gesdd/20603/19 "2019-02-10T18:07:36Z")

</div>

This is a simplified example taken out of context. I’m not assuming U==Vt. Most of the time my B will be complex valued and not symmetric/hermitian. But even in the simpler case I observe the qualitative difference, that’s why I simplified.

On a side note, I ended up using QR decompositions eventually, mainly for speed reasons.

---

<div class="post-metadata">

**Author:** ![Ralph\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ralph_smith/32/10344_2.png) [@Ralph\_Smith](https://discourse.julialang.org/u/Ralph_Smith)\
**Post date:** [February 10, 2019, 8:51pm UTC](https://discourse.julialang.org/t/svd-better-default-to-gesvd-instead-of-gesdd/20603/20 "2019-02-10T20:51:55Z")

</div>

On looking deeper, I see that LAPACK also has a Jacobi method for the SVD (as mentioned in [Calling Lapack's Jacobi SVD - dgesvj](https://discourse.julialang.org/t/calling-lapacks-jacobi-svd-dgesvj/20602)) This approach is not mentioned in the Users’ Guide, but (often?) gives singular values with _elementwise relative_ accuracy close to unit rounding error. It works for all of the s.v. for the companion matrices described above until N\>64 where things are so bad that `NaN`s are generated.

[Next page](https://discourse.julialang.org/t/svd-better-default-to-gesvd-instead-of-gesdd/20603.md?page=2)