# Truncated SVD of extended precision matrix

**URL:** <https://discourse.julialang.org/t/truncated-svd-of-extended-precision-matrix/78989>\
**Category:** Numerics\
**Tags:** extended-precision\
**Created:** [April 4, 2022, 8:03am UTC](https://discourse.julialang.org/t/truncated-svd-of-extended-precision-matrix/78989 "2022-04-04T08:03:22Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![Samuel3008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/samuel3008/32/25021_2.png) [@Samuel3008](https://discourse.julialang.org/u/Samuel3008)\
**Post date:** [April 4, 2022, 8:03am UTC](https://discourse.julialang.org/t/truncated-svd-of-extended-precision-matrix/78989/1 "2022-04-04T08:03:22Z")

</div>

Hi, is anybody aware of a way to compute the truncated SVD of an extended precision (128 bit floats in my case) matrix?  
For the floats, I’d like to use either MultiFloats.jl (seems faster in my testing, so this would be preferred) or DoubleFloats.jl. GenericLinearAlgebra.jl works to an extent but has two major problems that currently prevent me from using it: 1) It gets stuck on certain matrices ([GitHub issue](https://github.com/JuliaLinearAlgebra/GenericLinearAlgebra.jl/issues/81)) and 2) computes the entire, untruncated SVD, which is much more expensive.

I’m grateful for any hints or experiences you can share!

---

<div class="post-metadata">

**Author:** ![jd-foster](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jd-foster/32/35824_2.png) [@jd-foster](https://discourse.julialang.org/u/jd-foster)\
**Post date:** [April 5, 2022, 10:58am UTC](https://discourse.julialang.org/t/truncated-svd-of-extended-precision-matrix/78989/2 "2022-04-05T10:58:06Z")

</div>

I saw this earlier post:

> [@Truncated Singular Value Decomposition](https://discourse.julialang.org/t/truncated-singular-value-decomposition/42124):
>
> Hi. Given a matrix M I would like to compute its SVD truncated to rank k. I think this is possible without doing the full SVD. For example, Python has this: [sklearn.decomposition.TruncatedSVD — scikit-learn 1.1.2 documentation](https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.TruncatedSVD.html). The code I am currently using to do this is given below. The problem is that it computes SVD first, and then throws out the extra rows/columns, which can be quite costly if k is much smaller than rank of M. Is there a function in Julia to do this? using LinearAlgebra "re…

The recommendation is [TSVD.jl](https://github.com/JuliaLinearAlgebra/TSVD.jl).

---

<div class="post-metadata">

**Author:** ![Samuel3008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/samuel3008/32/25021_2.png) [@Samuel3008](https://discourse.julialang.org/u/Samuel3008)\
**Post date:** [April 5, 2022, 11:52am UTC](https://discourse.julialang.org/t/truncated-svd-of-extended-precision-matrix/78989/3 "2022-04-05T11:52:30Z")

</div>

Thanks, but unfortunately TSVD.jl calls LinearAlgebra‘s `svd`, and so doesn‘t work for types not supported by whatever BLAS julia uses.

---

<div class="post-metadata">

**Author:** ![andreasnoack](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andreasnoack/32/27_2.png) [@andreasnoack](https://discourse.julialang.org/u/andreasnoack)\
**Post date:** [April 5, 2022, 3:36pm UTC](https://discourse.julialang.org/t/truncated-svd-of-extended-precision-matrix/78989/4 "2022-04-05T15:36:11Z")

</div>

> [@Samuel3008](#):
>
> Thanks, but unfortunately TSVD.jl calls LinearAlgebra‘s `svd` , and so doesn‘t work for types not supported by whatever BLAS julia uses.

You should be able to load `GenericLinearAlgebra` which would extend `svd` to work with arbitrary element types.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [April 5, 2022, 4:40pm UTC](https://discourse.julialang.org/t/truncated-svd-of-extended-precision-matrix/78989/5 "2022-04-05T16:40:20Z")

</div>

> [@andreasnoack](#):
>
> You should be able to load `GenericLinearAlgebra` which would extend `svd` to work with arbitrary element types.

In particular, TSVD.jl calls [svd(::Bidiagonal)](https://github.com/JuliaLang/julia/blob/3fb132f0f514e98e54c962a03d2771e14cfdc947/stdlib/LinearAlgebra/src/bidiag.jl#L214-L221), and it looks like GenericLinearAlgebra.jl [provides a method for this](https://github.com/JuliaLinearAlgebra/GenericLinearAlgebra.jl/blob/bc746a2d7dc0c404376b7d4deca459bd2a21a331/src/svd.jl#L493-L499).

---

<div class="post-metadata">

**Author:** ![Samuel3008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/samuel3008/32/25021_2.png) [@Samuel3008](https://discourse.julialang.org/u/Samuel3008)\
**Post date:** [April 6, 2022, 11:50am UTC](https://discourse.julialang.org/t/truncated-svd-of-extended-precision-matrix/78989/6 "2022-04-06T11:50:32Z")

</div>

Thank you! I‘ll give this a try and report back.

---

<div class="post-metadata">

**Author:** ![Samuel3008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/samuel3008/32/25021_2.png) [@Samuel3008](https://discourse.julialang.org/u/Samuel3008)\
**Post date:** [April 30, 2022, 10:03am UTC](https://discourse.julialang.org/t/truncated-svd-of-extended-precision-matrix/78989/7 "2022-04-30T10:03:21Z")

</div>

Ok, so I finally got around to trying this:

```julia-repl
julia> using TSVD

julia> using GenericLinearAlgebra

julia> using DoubleFloats

julia> A = randn(Double64, 3, 5)
3×5 Matrix{Double64}:
  1.3801363771568418 1.4219282164212106 0.13499940802352295 -0.16125764163952927 1.6214322111600634
  0.04669843558325742 0.7916858087312778 -0.740636655635555 0.2954127906942305 -1.2628958884690384
 -0.9117030073895186 0.5428515183127985 -0.10123196960714188 0.9868982725916127 0.2833604981860745

julia> tsvd(A, 2)
ERROR: SVD for lower triangular bidiagonal matrices isn't implemented yet.
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:33
  [2] __svd!(B::LinearAlgebra.Bidiagonal{Double64, Vector{Double64}}, U::Matrix{Double64}, Vᴴ::Matrix{Double64}; tol::Double64, debug::Bool)
    @ GenericLinearAlgebra ~/.julia/packages/GenericLinearAlgebra/Bd8s0/src/svd.jl:196
  [3] _svd!(B::LinearAlgebra.Bidiagonal{Double64, Vector{Double64}}; tol::Double64, debug::Bool)
    @ GenericLinearAlgebra ~/.julia/packages/GenericLinearAlgebra/Bd8s0/src/svd.jl:232
  [4] #svd!#38
    @ ~/.julia/packages/GenericLinearAlgebra/Bd8s0/src/svd.jl:499 [inlined]
  [5] svd!
    @ ~/.julia/packages/GenericLinearAlgebra/Bd8s0/src/svd.jl:499 [inlined]
  [6] #svd#191
    @ /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/bidiag.jl:214 [inlined]
  [7] svd
    @ /Applications/Julia-1.7.app/Contents/Resources/julia/share/julia/stdlib/v1.7/LinearAlgebra/src/bidiag.jl:214 [inlined]
  [8] biLanczos(A::Matrix{Double64}, nvals::Int64; maxiter::Int64, initvec::Vector{Double64}, tolconv::Double64, tolreorth::Double64, stepsize::Int64, debug::Bool)
    @ TSVD ~/.julia/packages/TSVD/q42gG/src/svd.jl:175
  [9] _tsvd(A::Matrix{Double64}, nvals::Int64; maxiter::Int64, initvec::Vector{Double64}, tolconv::Double64, tolreorth::Double64, stepsize::Int64, debug::Bool)
    @ TSVD ~/.julia/packages/TSVD/q42gG/src/svd.jl:236
 [10] #tsvd#7
    @ ~/.julia/packages/TSVD/q42gG/src/svd.jl:345 [inlined]
 [11] tsvd(A::Matrix{Double64}, nvals::Int64)
    @ TSVD ~/.julia/packages/TSVD/q42gG/src/svd.jl:345
 [12] top-level scope
    @ REPL[5]:1

```

The relevant part in `GenericLinearAlgebra.jl` is this:

```julia
    if B.uplo === 'U'
        # actual code here, omitted for readability
    else
        # Just transpose the matrix
        error("SVD for lower triangular bidiagonal matrices isn't implemented yet.")
    end

```

Seems like this would be easiest to fix in `TSVD.jl` by, as the comment suggests, transposing the bidiagonal matrix before feeding it to `GenericLinearAlgebra.jl`, right?