Inverse with destination?


#1

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?


#2

Are you looking for

dest = inv(X)

?


#3

You may be able to use

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.


#4

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.


#5

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.


#6

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.


#7

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.


#8

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


#9

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: 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:

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⁻¹