# Transformations in ArnoldiMethod.jl
**URL:** https://discourse.julialang.org/t/transformations-in-arnoldimethod-jl/109953
**Category:** General Usage
**Created:** [February 8, 2024, 9:48pm UTC](https://discourse.julialang.org/t/transformations-in-arnoldimethod-jl/109953 "2024-02-08T21:48:28Z")
**Posts on this page:** 3
**Page:** 1
### Author:  [@densmojd](https://discourse.julialang.org/u/densmojd)
#### Post date: [February 8, 2024, 9:48pm UTC](https://discourse.julialang.org/t/transformations-in-arnoldimethod-jl/109953/1 "2024-02-08T21:48:28Z")
I have a question about the documentation of “transformations” (i.e., multiplying the inverse of a matrix times a vector without explicitly forming the inverse) in [ArnoldiMethod.jl](https://github.com/JuliaLinearAlgebra/ArnoldiMethod.jl), although I think the topic is of broader interest than just this package.
From the [documentation](https://julialinearalgebra.github.io/ArnoldiMethod.jl/stable/usage/02_spectral_transformations.html):
```julia
using ArnoldiMethod, LinearAlgebra, LinearMaps
# Define a matrix whose eigenvalues you want
A = rand(100,100)
# Factorizes A and builds a linear map that applies inv(A) to a vector.
function construct_linear_map(A)
F = factorize(A)
LinearMap{eltype(A)}((y, x) -> ldiv!(y, F, x), size(A,1), ismutating=true)
end
# Target the largest eigenvalues of the inverted problem
decomp, = partialschur(construct_linear_map(A), nev=4, tol=1e-5, restarts=100, which=LM())
λs_inv, X = partialeigen(decomp)
```
However, I would replace the definition of `construct_linear_map` and the call to `partial_schur` with the somewhat simpler
```julia
F = factorize(A)
linearmap = LinearMap{eltype(A)}((y, x) -> ldiv!(y, F, x), size(A,1), ismutating=true)
decomp, = partialschur(linearmap, nev=4, tol=1e-5, restarts=100, which=LM())
```
The example for a generalized eigenvalue problem is more complicated:
```julia
using ArnoldiMethod, LinearAlgebra, LinearMaps
# Define the matrices of the generalized eigenvalue problem
A, B = rand(100,100), rand(100,100)
struct ShiftAndInvert{TA,TB,TT}
A_lu::TA
B::TB
temp::TT
end
function (M::ShiftAndInvert)(y,x)
mul!(M.temp, M.B, x)
ldiv!(y, M.A_lu, M.temp)
end
function construct_linear_map(A,B)
a = ShiftAndInvert(factorize(A),B,Vector{eltype(A)}(undef, size(A,1)))
LinearMap{eltype(A)}(a, size(A,1), ismutating=true)
end
# Target the largest eigenvalues of the inverted problem
decomp, = partialschur(construct_linear_map(A, B), nev=4, tol=1e-5, restarts=100, which=LM())
λs_inv, X = partialeigen(decomp)
```
I would simplify this code by replacing the definitions of `ShiftAndInvert` and `construct_linear_map` and the call to `partialschur` with
```julia
F = factorize(A)
linearmap = LinearMap{eltype(A)}((y, x) -> ldiv!(y, F, B * x), size(A,1), ismutating=true)
decomp, = partialschur(linearmap, nev=4, tol=1e-5, restarts=100, which=LM())
```
When I use my approach in my own work (full disclosure: not these eigenvalue problems specifically), it seems to work fine (e.g., `@code_warntype` doesn’t report any type instabilities). So why what appears to me more complicated code?
---
### Author:  [@abraemer](https://discourse.julialang.org/u/abraemer)
#### Post date: [February 9, 2024, 8:02am UTC](https://discourse.julialang.org/t/transformations-in-arnoldimethod-jl/109953/2 "2024-02-09T08:02:57Z")
Welcome to the discourse 🙂
In the first example, I don’t understand why you say that your version is simpler. To me the code does the very same thing, except that defining a function to construct the LinearMap has two advantages:
1. You don’t need to recall all the things you put into the LinearMap everytime you want to find some eigenvalues of a matrix
2. Your solution should be a bit slower if it lives in global scope since you repeatedly access `F` and [that is slow](https://docs.julialang.org/en/v1/manual/performance-tips/#Avoid-untyped-global-variables). If it lives inside a larger function and you override `F` later, you could still run into [performance issues due to captured variables](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-captured). Which is much more unlikely if it is wrapped up like in the example code.
Your second solution has the same “maybe-problems” with repeated access to `F` and `B` and on top of that causes allocations on every application of the LinearMap because `B*x` allocates a new vector. The example avoids that by preallocating some memory and using it to store the result of the multiplication with `mul!`.
You are of course right, that there are different way of achieving the same thing and what is “easiest” or “more complicated” might not be the same for everyone. To me both codes read quite clear. The example code just tied up the logic into a type and neat construction function for the LinearMap. E.g. you can have the same thing without a struct:
```julia
linearmap = let F = factorize(A),
B = B, # to avoid performance issues from captured variables
temp = Vector{eltype(A)}(undef, size(A,1)
LinearMap{eltype(A)}((y, x) -> ldiv!(y, F, mul!(temp, B, x)), size(A,1), ismutating=true)
end
decomp, = partialschur(linearmap, nev=4, tol=1e-5, restarts=100, which=LM())
```
---
### Author:  [@densmojd](https://discourse.julialang.org/u/densmojd)
#### Post date: [February 9, 2024, 7:21pm UTC](https://discourse.julialang.org/t/transformations-in-arnoldimethod-jl/109953/3 "2024-02-09T19:21:55Z")
Thanks for the welcome!
All your points are good. For what it’s worth, I thought there might be a global scope/capture variable issue, but in my own code all of this is buried in a function and `@code_warntype` seemed happy. What I didn’t realize, though, that `B*x` allocates; thanks for pointing that out.