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

Internals & Design
17

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
product_chain_stable2000×700 78.5 KB

MWE:

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)