# In-place svd?

**URL:** <https://discourse.julialang.org/t/in-place-svd/60567>\
**Category:** Performance\
**Created:** [May 5, 2021, 11:52am UTC](https://discourse.julialang.org/t/in-place-svd/60567 "2021-05-05T11:52:12Z")\
**Posts on this page:** 17\
**Page:** 1

<div class="post-metadata">

**Author:** ![tomtom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tomtom/32/5106_2.png) [@tomtom](https://discourse.julialang.org/u/tomtom)\
**Post date:** [May 5, 2021, 11:52am UTC](https://discourse.julialang.org/t/in-place-svd/60567/1 "2021-05-05T11:52:12Z")

</div>

I need to do SVD many many times, so I wanna have a in-place version of svd().

What I mean by “in-place” is like:

```julia
svd!(U, S, V, A)

```

which works out the svd of `A` and _ **writes the output matrices into the pre-allocated `U`, `S` and `V`** _.

Note the `LinearAlgebra.svd!`, `LinearAlgebra.LAPACK.ggsvd!` and `LinearAlgebra.LAPACK.gesvd!` are all **reusing** the memory of input for immediate calculations, but **all outputs are allocated**.

Thanks.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [May 5, 2021, 1:17pm UTC](https://discourse.julialang.org/t/in-place-svd/60567/2 "2021-05-05T13:17:44Z")

</div>

In general, this won’t help that much for big matrices since svd is O(n^3), and for small matrices, you probably want to be using StaticArrays. That said, it wouldn’t hurt to add this.

---

<div class="post-metadata">

**Author:** ![tomtom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tomtom/32/5106_2.png) [@tomtom](https://discourse.julialang.org/u/tomtom)\
**Post date:** [May 5, 2021, 1:52pm UTC](https://discourse.julialang.org/t/in-place-svd/60567/3 "2021-05-05T13:52:03Z")

</div>

Google a bit found that recently a paper claims O(n^2) for a deterministic SVD algorithm.

My point is, the `!()` functions in `LinearAlgebra` (e.g. `eigen!()`) are not “in-place” in conventional sense. And “real in-place” functions are very needed!

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [May 5, 2021, 1:56pm UTC](https://discourse.julialang.org/t/in-place-svd/60567/4 "2021-05-05T13:56:05Z")

</div>

I believe that the O(n^2) is only for very specific structures of matrices (eg banded with fixed band size). Real in-place methods wouldn’t be bad, but won’t lead to major speedups.

---

<div class="post-metadata">

**Author:** ![jling](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jling/32/212909_2.png) [@jling](https://discourse.julialang.org/u/jling)\
**Post date:** [May 5, 2021, 1:56pm UTC](https://discourse.julialang.org/t/in-place-svd/60567/5 "2021-05-05T13:56:43Z")

</div>

[docs](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.svd!)

am I missing something? there seems to be a `svd!(A)` method that is just what you want?

---

<div class="post-metadata">

**Author:** ![juthohaegeman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juthohaegeman/32/8620_2.png) [@juthohaegeman](https://discourse.julialang.org/u/juthohaegeman)\
**Post date:** [May 5, 2021, 2:17pm UTC](https://discourse.julialang.org/t/in-place-svd/60567/6 "2021-05-05T14:17:00Z")

</div>

By studying the Julia implementation of `LAPACK.gesvd!` or `LAPACK.gesdd!` and the LAPACK manual, you can easily write such a function `svd!(U, S, V, A)` that does exactly what you want. You’ll probably even want to provide an additional argument for a work buffer that you can also recycle in between iterations.

Note that there is no option in LAPACK not to destroy the contents of `A`. Depending on the `jobu` and `jobv` parameters, you can recycle the memory of `A` to store the final `U` or `V`, so that you do not need to allocate one of those, but even if you have allocated space for those, you cannot keep the contents of `A`. So you would need to take a copy of `A` if you still need to use it afterwards.

---

<div class="post-metadata">

**Author:** ![tomtom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tomtom/32/5106_2.png) [@tomtom](https://discourse.julialang.org/u/tomtom)\
**Post date:** [May 5, 2021, 3:25pm UTC](https://discourse.julialang.org/t/in-place-svd/60567/7 "2021-05-05T15:25:37Z")

</div>

I guess speed up would be significant, especially if lots of garbage collections are needed for running non in-place versions many many times.

---

<div class="post-metadata">

**Author:** ![tomtom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tomtom/32/5106_2.png) [@tomtom](https://discourse.julialang.org/u/tomtom)\
**Post date:** [May 5, 2021, 3:28pm UTC](https://discourse.julialang.org/t/in-place-svd/60567/8 "2021-05-05T15:28:01Z")

</div>

`LinearAlgebra.svd!()` allocates memory for the output, i.e., I could NOT pre-allocate the memory for output. It just reuse the memory of input for immediate calculations.

---

<div class="post-metadata">

**Author:** ![jling](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jling/32/212909_2.png) [@jling](https://discourse.julialang.org/u/jling)\
**Post date:** [May 5, 2021, 3:33pm UTC](https://discourse.julialang.org/t/in-place-svd/60567/9 "2021-05-05T15:33:36Z")

</div>

> [@tomtom](#):
>
> which works out the svd of `A` and _ **writes the output matrices into the pre-allocated `U` , `S` and `V`** _ .

> [@tomtom](#):
>
> I could NOT pre-allocate the memory for output.

I’m confused to what you want exactly

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [May 5, 2021, 3:38pm UTC](https://discourse.julialang.org/t/in-place-svd/60567/10 "2021-05-05T15:38:50Z")

</div>

The reason this isn’t significant is that `svd` is slow. Specifically, for a 100x100 matrix, `svd` takes 1.5ms, while allocating arrays for the output is about 7us, and a gc run is in the 10ms range, but won’t be called frequently. As such, for a loop of running `svd` on a 100x100 matrix, I’m seeing .3% gc time on average. For smaller matrices, it can matter more, but I’m only seeing 1% gc time for 20x20 matrices, and if you get much smaller, you should be using `StaticArrays` anyway.

---

<div class="post-metadata">

**Author:** ![tomtom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tomtom/32/5106_2.png) [@tomtom](https://discourse.julialang.org/u/tomtom)\
**Post date:** [May 6, 2021, 11:57am UTC](https://discourse.julialang.org/t/in-place-svd/60567/11 "2021-05-06T11:57:50Z")

</div>

what about the complexity of `LinearAlgebra.eigen()` then? is it worth to develop in-place version of it?  
Thanks.

---

<div class="post-metadata">

**Author:** ![juthohaegeman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juthohaegeman/32/8620_2.png) [@juthohaegeman](https://discourse.julialang.org/u/juthohaegeman)\
**Post date:** [May 6, 2021, 12:08pm UTC](https://discourse.julialang.org/t/in-place-svd/60567/12 "2021-05-06T12:08:39Z")

</div>

`eigen` is slower than `svd`, unless its a hermitian matrix. Probably @Oscar_Smith is correct for factorisations. For multiplying matrices, even though the multiplication is also O(N^3), while the required memory for the result is O(N^2), I do think it can make a difference, just because this O(N^3) has been so optimised (and has a much smaller prefactor than in the svd or eigen case).

It has been my experience that the Julia’s GC, despite certainly having been much improved, still has some overhead when a lot of big temporary arrays are involved, in comparison to Matlab or Python. But I would need to do more thorough benchmarking and find a good test case; it’s hard to carefully compare this in a way that excludes all other possible causes. It is however a fact that in TensorOperations.jl I added a package wide cache structure specifically for storing intermediate arrays in tensor contractions, and (provided they all fit in memory), this can improve run time with a significant fraction.

---

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [May 24, 2022, 7:59am UTC](https://discourse.julialang.org/t/in-place-svd/60567/14 "2022-05-24T07:59:13Z")

</div>

The use of `!` in Julia is a little confusing.  
For low level programmers functions which overwrite the input (In place) are also usually non allocating while in Julia that’s not the case.

It would be great to have an option for non allocating API for BLAS / LAPACK operations.  
Especially in the context of advancement in `StaticCompiler.jl`.  
If one wants to build a static / dynamic library using Julia, it would be great to be able to make sure it can be used with no allocations at all.

---

<div class="post-metadata">

**Author:** ![moeddel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/moeddel/32/18641_2.png) [@moeddel](https://discourse.julialang.org/u/moeddel)\
**Post date:** [May 24, 2022, 12:18pm UTC](https://discourse.julialang.org/t/in-place-svd/60567/15 "2022-05-24T12:18:41Z")

</div>

If I recall correctly, the `!` does not make any statements about allocations. It is merely a hint that [a function modifies its arguments](https://docs.julialang.org/en/v1/manual/style-guide/#bang-convention).

---

<div class="post-metadata">

**Author:** ![RoyiAvital](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/royiavital/32/571_2.png) [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)\
**Post date:** [May 24, 2022, 12:25pm UTC](https://discourse.julialang.org/t/in-place-svd/60567/16 "2022-05-24T12:25:21Z")

</div>

Indeed, hence I said for people coming form low level programming. When you use functions which operate on the given buffers they don’t allocate.

Hence I wish for such API’s in Julia as well. My dream is to have an option to solve Linear System with zero allocations. Both for Sparse and Dense matrices.

---

<div class="post-metadata">

**Author:** ![moeddel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/moeddel/32/18641_2.png) [@moeddel](https://discourse.julialang.org/u/moeddel)\
**Post date:** [May 24, 2022, 7:24pm UTC](https://discourse.julialang.org/t/in-place-svd/60567/17 "2022-05-24T19:24:45Z")

</div>

Agreed that it would be nice to have.

On the other hand you could also code this on your own. The implementation of [gesvd](https://github.com/JuliaLang/julia/blob/742b9abb4dd4621b667ec5bb3434b8b3602f96fd/stdlib/LinearAlgebra/src/lapack.jl#L1560) is quite straight forward and you can easily see, where the allocations happen.

---

<div class="post-metadata">

**Author:** ![Algopaul](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/algopaul/32/37603_2.png) [@Algopaul](https://discourse.julialang.org/u/Algopaul)\
**Post date:** [July 2, 2022, 1:19pm UTC](https://discourse.julialang.org/t/in-place-svd/60567/18 "2022-07-02T13:19:07Z")

</div>

If you still have this problem, you may want to check out [KWLinalg](https://github.com/Algopaul/KWLinalg)  
This is a package, in which (among other things) functors are provided that already include the necessary workspace and output matrices. Say you have a matrix

```julia
A = rand(ComplexF64, 5, 4)

```

Then you can define the functor

```julia
f = svd_functor_divconquer(5, 4, ComplexF64)

```

After that, you can call

```julia
U, S, Vt = f(A)

```

with no further allocations (unless you are in the main-scope). See the example in the README of [KWLinalg](https://github.com/Algopaul/KWLinalg) for a benchmark-test.