# Inverse with destination?

**URL:** <https://discourse.julialang.org/t/inverse-with-destination/2267>\
**Category:** General Usage\
**Created:** [February 24, 2017, 7:51am UTC](https://discourse.julialang.org/t/inverse-with-destination/2267 "2017-02-24T07:51:59Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Andrea\_Pagnani](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andrea_pagnani/32/3136_2.png) [@Andrea\_Pagnani](https://discourse.julialang.org/u/Andrea_Pagnani)\
**Post date:** [February 24, 2017, 7:51am UTC](https://discourse.julialang.org/t/inverse-with-destination/2267/1 "2017-02-24T07:51:59Z")

</div>

Dear all,

I might be wrong, but in the plethora of BLAS and LAPACK functions, I was not able to find a function of the type

> inv(dest, X)

where X is the matrix to invert and dest is the destination.  
Any suggestion?

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [February 24, 2017, 8:12am UTC](https://discourse.julialang.org/t/inverse-with-destination/2267/2 "2017-02-24T08:12:16Z")

</div>

Are you looking for

dest = inv(X)

?

---

<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 24, 2017, 9:52am UTC](https://discourse.julialang.org/t/inverse-with-destination/2267/3 "2017-02-24T09:52:35Z")

</div>

You may be able to use

```julia
Base.LinAlg.LAPACK.getri!

```

but note that it overwrites the argument, so you may want to copy that first. I am not aware of a LAPACK routine that puts the inverse or the factorization in another matrix, they usually overwrite when that makes sense.

---

<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:** [February 24, 2017, 12:39pm UTC](https://discourse.julialang.org/t/inverse-with-destination/2267/4 "2017-02-24T12:39:14Z")

</div>

You can use `A_ldiv_B!(A, X)` after initializing the contents of `X` to the identity matrix.

However, I would first make sure that you really want to compute the inverse matrix. Most of the time, you should just compute `LU = lufact(A)` (or `lufact!` to do it in-place) — almost anything you might want to do with the inverse matrix can be done about as efficiently with the LU factorization, which saves you the added expense of computing the inverse.

---

<div class="post-metadata">

**Author:** ![Andrea\_Pagnani](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andrea_pagnani/32/3136_2.png) [@Andrea\_Pagnani](https://discourse.julialang.org/u/Andrea_Pagnani)\
**Post date:** [February 25, 2017, 10:23am UTC](https://discourse.julialang.org/t/inverse-with-destination/2267/5 "2017-02-25T10:23:29Z")

</div>

> julia\> N=3; dest = zeros(N,N); p1 = pointer(dest); dest=inv(rand(N,N)); p2 = pointer(dest); p1==p2  
> false

It is not what I’m looking for as it creates a copy.

---

<div class="post-metadata">

**Author:** ![Andrea\_Pagnani](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andrea_pagnani/32/3136_2.png) [@Andrea\_Pagnani](https://discourse.julialang.org/u/Andrea_Pagnani)\
**Post date:** [February 25, 2017, 2:16pm UTC](https://discourse.julialang.org/t/inverse-with-destination/2267/6 "2017-02-25T14:16:07Z")

</div>

As the matrix to be inverted is posdef the following (thanks to @carlobaldassi)

> function inplaceinverse!(dest,source)  
> copy!(dest, source)  
> Base.LinAlg.inv!(cholfact!(dest))  
> end

suits my needs. Thanks to @stevengj for his suggestion but unfortunately I need the diagonal of the inverse of a posdef matrix. Maybe, but this is a question on linear algebra, there is a smarter strategy to obtain the diagonal of the inverse of a posdef matrix in a faster and less memory requiring way.

---

<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 25, 2017, 2:57pm UTC](https://discourse.julialang.org/t/inverse-with-destination/2267/7 "2017-02-25T14:57:22Z")

</div>

> [@Andrea\_Pagnani](#):
>
> Maybe, but this is a question on linear algebra, there is a smarter strategy to obtain the diagonal of the inverse of a posdef matrix in a faster and less memory requiring way.

That sometimes happens when the positive definite matrix was obtained from some other matrix, and you work with the latter directly, eg using a factorization.

---

<div class="post-metadata">

**Author:** ![StefanKarpinski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stefankarpinski/32/24_2.png) [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)\
**Post date:** [March 2, 2017, 4:12pm UTC](https://discourse.julialang.org/t/inverse-with-destination/2267/8 "2017-03-02T16:12:52Z")

</div>

Would it make sense to have an `inv!(dest, X)` function in Base?

---

<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:** [March 2, 2017, 10:08pm UTC](https://discourse.julialang.org/t/inverse-with-destination/2267/9 "2017-03-02T22:08:44Z")

</div>

> [@Andrea\_Pagnani](#):
>
> Maybe, but this is a question on linear algebra, there is a smarter strategy to obtain the diagonal of the inverse of a posdef matrix in a faster and less memory requiring way.

There are iterative methods for this, but they are most beneficial when your matrix is sparse (or has some other special structure that lets you compute matrix × vector quickly).

See also this discussion: [matrices - Find diagonal of inverse matrix - Mathematics Stack Exchange](http://math.stackexchange.com/questions/978051/find-diagonal-of-inverse-matrix) … with direct methods for dense matrices like Cholesky, it doesn’t seem like you can save more than a factor of two compared to computing the whole inverse matrix, and in practice it will probably be less than that because of the difficulty of beating BLAS.

For an m×m SPD matrix `A`, pre-allocated matrices `Acopy` and `X` of the same size, and a preallocated array `d` of length `m`, a straightforward and fast way to extract the diagonal of the inverse in-place is:

```julia
Acopy .= A # make a copy if you don't want to overwrite A
X .= 0
for i = 1:m; X[i,i] = 1; end # set X to identity matrix
C = cholfact!(Acopy)
A_ldiv_B!(C, X)
for i = 1:m; d[i] = X[i,i]; end # extract diagonal of A⁻¹

```

---

<div class="post-metadata">

**Author:** ![e3c6](https://avatars.discourse-cdn.com/v4/letter/e/e79b87/32.png) [@e3c6](https://discourse.julialang.org/u/e3c6)\
**Post date:** [May 17, 2019, 4:43pm UTC](https://discourse.julialang.org/t/inverse-with-destination/2267/10 "2019-05-17T16:43:03Z")

</div>

> [@StefanKarpinski](#):
>
> Would it make sense to have an `inv!(dest, X)` function in Base?

I think it would be useful to have `inv!(dest, X)`. I don’t know if after two years there is now another way to do this?

---

<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:** [May 17, 2019, 6:32pm UTC](https://discourse.julialang.org/t/inverse-with-destination/2267/11 "2019-05-17T18:32:25Z")

</div>

> [@e3c6](#):
>
> I don’t know if after two years there is now another way to do this?

Similar to [this post](https://discourse.julialang.org/t/inverse-with-destination/2267/6), there is `LinearAlgebra.inv!` (currently unexported and undocumented), which works on factorization objects. So you can do:

```julia
LinearAlgebra.inv!(lu!(A))

```

for example, to perform in-place inversion. Or use `cholesky!` instead of `lu!` if you have an SPD matrix. If you want a different destination, use `copyto!(dest, A)` instead of `A`.

---

<div class="post-metadata">

**Author:** ![e3c6](https://avatars.discourse-cdn.com/v4/letter/e/e79b87/32.png) [@e3c6](https://discourse.julialang.org/u/e3c6)\
**Post date:** [May 17, 2019, 9:28pm UTC](https://discourse.julialang.org/t/inverse-with-destination/2267/12 "2019-05-17T21:28:41Z")

</div>

I’m computing the inverse of a `Tridiagonal` matrix, and I want to place it in a pre-allocated array. In this case I cannot use in-place inverse. I would need `inv!(dest, X)` or something equivalent. Any suggestions?

---

<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:** [May 18, 2019, 2:45am UTC](https://discourse.julialang.org/t/inverse-with-destination/2267/13 "2019-05-18T02:45:31Z")

</div>

> [@e3c6](#):
>
> I’m computing the inverse of a `Tridiagonal` matrix

Don’t do this. Use the LU factorization for whatever you are trying to do.

You can solve NxN tridiagonal systems in O(N) operations. The inverse, on the other hand, is a dense matrix that throws this advantage away.

---

<div class="post-metadata">

**Author:** ![e3c6](https://avatars.discourse-cdn.com/v4/letter/e/e79b87/32.png) [@e3c6](https://discourse.julialang.org/u/e3c6)\
**Post date:** [May 18, 2019, 8:19am UTC](https://discourse.julialang.org/t/inverse-with-destination/2267/14 "2019-05-18T08:19:45Z")

</div>

> [@stevengj](#):
>
> Don’t do this.

I am doing a weighted linear least squares regression and need to compute error bars. Unfortunately at some point there is a tridiagonal matrix that needs to be inverted (and I don’t see a way around this inversion).

---

<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:** [May 18, 2019, 8:34am UTC](https://discourse.julialang.org/t/inverse-with-destination/2267/15 "2019-05-18T08:34:36Z")

</div>

> [@e3c6](#):
>
> I don’t see a way around this inversion

Unless you are using A^{-1} directly, any A^{-1}b or similar is best implemented as `A \ b`.

Calculating an inverse of a matrix is the numerical equivalent of code smell emanating from a medium-size trout that has been under the sofa for two weeks at 20°C.

---

<div class="post-metadata">

**Author:** ![e3c6](https://avatars.discourse-cdn.com/v4/letter/e/e79b87/32.png) [@e3c6](https://discourse.julialang.org/u/e3c6)\
**Post date:** [May 18, 2019, 8:50am UTC](https://discourse.julialang.org/t/inverse-with-destination/2267/16 "2019-05-18T08:50:40Z")

</div>

> [@Tamas\_Papp](#):
>
> Unless you are using A^{-1} directly

I am.

If you are interested, this is what I’m trying to do: [Weighted least squares - Wikipedia](https://en.wikipedia.org/wiki/Weighted_least_squares#Parameter_errors_and_correlation)

---

<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:** [May 18, 2019, 8:59am UTC](https://discourse.julialang.org/t/inverse-with-destination/2267/17 "2019-05-18T08:59:02Z")

</div>

Never use OLS formulas for computation. For WLS, there are much better options, eg premultiplying by the relevant Cholesky factor or its inverse (depending on details and how you obtain W).

---

<div class="post-metadata">

**Author:** ![e3c6](https://avatars.discourse-cdn.com/v4/letter/e/e79b87/32.png) [@e3c6](https://discourse.julialang.org/u/e3c6)\
**Post date:** [May 18, 2019, 9:04am UTC](https://discourse.julialang.org/t/inverse-with-destination/2267/18 "2019-05-18T09:04:14Z")

</div>

> [@Tamas\_Papp](#):
>
> (depending on details and how you obtain W)

W is the inverse of the covariance matrix of observation errors (errors are not independent). How you suggest I compute the covariance matrix of the parameters (denoted M^\beta in the Wikipedia link)? (I only need the variances).

---

<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:** [May 18, 2019, 2:01pm UTC](https://discourse.julialang.org/t/inverse-with-destination/2267/19 "2019-05-18T14:01:25Z")

</div>

> [@e3c6](#):
>
> How you suggest I compute the covariance matrix of the parameters (denoted M^\beta in the Wikipedia link)?

Don’t — computer the Cholesky factor of M^\beta. If you have W=LL', it is fairly simple matrix algebra to work it out, but it may be better to start with the QR decomposition of the demeaned observations.

---

<div class="post-metadata">

**Author:** ![e3c6](https://avatars.discourse-cdn.com/v4/letter/e/e79b87/32.png) [@e3c6](https://discourse.julialang.org/u/e3c6)\
**Post date:** [May 19, 2019, 9:13pm UTC](https://discourse.julialang.org/t/inverse-with-destination/2267/20 "2019-05-19T21:13:47Z")

</div>

I have read many times the recommendation to avoid matrix inversion and use `LU` instead (or whatever factorization is applicable). However in the particular case of linear regression, it is always pointed out that if parameter errors are variances, the inversion is generally unavoidable.

For example, from [Wikipedia](https://en.wikipedia.org/wiki/Invertible_matrix#Applications) (there might be better sources but I’m too lazy to find them now):

> Although an explicit inverse is not necessary to estimate the vector of unknowns, it is unavoidable to estimate their precision, found in the diagonal of the posterior covariance matrix of the vector of unknowns.

Also in general parameter estimation by maximum likelihood, the inverse of the Hessian of the likelihood allows one to estimate parameter errors.

As far as I am aware, these are applications where matrix inversion is unavoidable (but please correct me if I’m wrong).

I just want to make the case that `inv!(dest, A)` could be useful in situations like this (specially if `A` is sparse, so in-place `inv!(A)` is not useful), and I suspect it should not be so hard to have?

[Next page](https://discourse.julialang.org/t/inverse-with-destination/2267.md?page=2)