# Matrix (linear map function) exponential multiplied with vector in ExponentialUtilities.jl

**URL:** <https://discourse.julialang.org/t/matrix-linear-map-function-exponential-multiplied-with-vector-in-exponentialutilities-jl/126957>\
**Category:** General Usage\
**Tags:** linearalgebra, physics, sparse, matrix\
**Created:** [March 14, 2025, 3:21pm UTC](https://discourse.julialang.org/t/matrix-linear-map-function-exponential-multiplied-with-vector-in-exponentialutilities-jl/126957 "2025-03-14T15:21:38Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![eveningsilverfox](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eveningsilverfox/32/51932_2.png) [@eveningsilverfox](https://discourse.julialang.org/u/eveningsilverfox)\
**Post date:** [March 14, 2025, 3:21pm UTC](https://discourse.julialang.org/t/matrix-linear-map-function-exponential-multiplied-with-vector-in-exponentialutilities-jl/126957/1 "2025-03-14T15:21:38Z")

</div>

I am trying to multiply the exponential of a matrix (available as a function which linearly maps the input vector to an output vector), to a given vector, using expv from ExponentialUtilities.jl

```julia
using MKL
using LinearAlgebra
using ExponentialUtilities

function matmap(x,aa,bb)
    Nx = size(x,1);
    y = zeros(eltype(x), Nx);
    for gh = 1:Nx
        @views y[gh] += aa*x[gh];
        if gh>1 && gh<Nx
            @views y[gh] += -bb*x[gh-1]-bb*x[gh+1];
        end
        @views y[1] += -bb*x[2];
        @views y[Nx] += -bb*x[Nx-1];
    end

    return y
end

dt = 0.01;
aa = 0.9; bb = 0.5;
Nx = 1000;
xin = rand(ComplexF64, Nx);
expv(-im*dt,x->matmap(x,aa,bb),xin);

```

I am getting the error

```julia
MethodError: no method matching size(::var"#7#8", ::Int64)
The function `size` exists, but no method is defined for this combination of argument types.

```

Any ideas on how to implement it correctly? In reality I deal with massive matrices, which I must apply with the function above as they are very sparse. So I cannot use the actual matrices. The example provided here is only a representative case, not the actual matrix I use.

---

<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 14, 2025, 5:27pm UTC](https://discourse.julialang.org/t/matrix-linear-map-function-exponential-multiplied-with-vector-in-exponentialutilities-jl/126957/2 "2025-03-14T17:27:45Z")

</div>

> [@eveningsilverfox](#):
>
> Any ideas on how to implement it correctly? In reality I deal with massive matrices, which I must apply with the function above as they are very sparse. So I cannot use the actual matrices. The example provided here is only a representative case, not the actual matrix I use.

I think it’s expecting a “matrix-like” object that supports a `size` method, not just a function. However, you can wrap a “matrix-like” object around your function with LinearMaps.jl, in particular with `LinearMap(x->matmap(x,aa,bb), Nx)` which creates an object that has `size == (Nx, Nx)`.

---

<div class="post-metadata">

**Author:** ![eveningsilverfox](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eveningsilverfox/32/51932_2.png) [@eveningsilverfox](https://discourse.julialang.org/u/eveningsilverfox)\
**Post date:** [March 14, 2025, 9:50pm UTC](https://discourse.julialang.org/t/matrix-linear-map-function-exponential-multiplied-with-vector-in-exponentialutilities-jl/126957/3 "2025-03-14T21:50:37Z")

</div>

I think LinearMap makes things slow when the function creating the map takes as arguments very large arrays (in my physics simulation, the Hamiltonian matrix/map maps basis vectors in a very large space, whose labels require hashing and recalculation using the stored basis states in the large array). My allocations are through the roof.  
Creating the sparse matrix once and feeding it makes it run faster.

Is there a way to avoid massive allocations with LinearMap in this case?

---

<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 14, 2025, 9:55pm UTC](https://discourse.julialang.org/t/matrix-linear-map-function-exponential-multiplied-with-vector-in-exponentialutilities-jl/126957/4 "2025-03-14T21:55:11Z")

</div>

> [@eveningsilverfox](#):
>
> Is there a way to avoid massive allocations with LinearMap in this case?

It sounds like the allocations have nothing to do with `LinearMap` and everything to do with the implementation of your function `matmap`? It’s hard to say how to make that function faster without seeing what it looks like, since it sounds like you are talking about a different function than the example you gave above?

(Note that you can reduce allocations in `LinearMap` by providing a mutating function `f!(y, x)` that writes the result of the map in-place into `y`, and then passing `ismutating=true` to the `LinearMap` constructor. In this case the allocations should be purely a question of how you implemented `f!`.)

---

<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 14, 2025, 10:10pm UTC](https://discourse.julialang.org/t/matrix-linear-map-function-exponential-multiplied-with-vector-in-exponentialutilities-jl/126957/6 "2025-03-14T22:10:05Z")

</div>

> [@eveningsilverfox](#):
>
> `ind1 = partialsortperm(abs.(tagssort .- tag1), 1);`

Can’t you cache your `ind2` values somehow rather than doing this for every application of the matrix? You probably need to think about your data structures more carefully here.

---

<div class="post-metadata">

**Author:** ![eveningsilverfox](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eveningsilverfox/32/51932_2.png) [@eveningsilverfox](https://discourse.julialang.org/u/eveningsilverfox)\
**Post date:** [March 14, 2025, 10:12pm UTC](https://discourse.julialang.org/t/matrix-linear-map-function-exponential-multiplied-with-vector-in-exponentialutilities-jl/126957/7 "2025-03-14T22:12:06Z")

</div>

You are correct, since the matrix/map follows a definite tridiagonal-like pattern, the indices `ind1` and `ind2`can be precomputed and passed in another array.

But this element is common for both the matrix-based `Hammat` and map-based `Hammap`…

---

<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 14, 2025, 10:15pm UTC](https://discourse.julialang.org/t/matrix-linear-map-function-exponential-multiplied-with-vector-in-exponentialutilities-jl/126957/8 "2025-03-14T22:15:12Z")

</div>

> [@eveningsilverfox](#):
>
> You are correct, since the matrix/map follows a definite tridiagonal-like pattern,

If the matrix is that structured, I’m confused about why you need to store indices (or compute them with complicated data structures) at all. The whole point of “matrix-free” algorithms (where you _don’t_ allocate an explicit sparse-matrix data structure with arrays of indices) is to exploit this kind of regularity.

(By the way, `partialsortperm(array, 1)` is equivalent to `findmin`, no?)

---

<div class="post-metadata">

**Author:** ![eveningsilverfox](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eveningsilverfox/32/51932_2.png) [@eveningsilverfox](https://discourse.julialang.org/u/eveningsilverfox)\
**Post date:** [March 14, 2025, 10:18pm UTC](https://discourse.julialang.org/t/matrix-linear-map-function-exponential-multiplied-with-vector-in-exponentialutilities-jl/126957/9 "2025-03-14T22:18:55Z")

</div>

By “tridiagonal-like”, I mean it corresponds to nearest neighbour connections on a lattice. But the basis states correspond to an arbitrary number of bosons sitting on each lattice site. So in the ordered basis, the Hamiltonian is not tridiagonal at all as whether any two basis elements actually have a non-zero matrix element between them depends on the number of bosons on each site given by the basis elements. So while it is definite, I can not hardcode the Hamiltonian matrix.  
Nevertheless, since we visit these elements in exactly the same way, I can store the ind1,ind2 from a single pass in a big vector, and reuse them.

Yes, I should use findmin instead. The partialsortperm is a vestigial code from many other stuff I have been copying from.

Also, when I benchmark just the Hamiltonians, the map is actually faster. It is only slower when I feed it to partialschur using a LinearMap.

```julia
instate = rand(Float64, Nstate);
@btime Hammat($zeta,$U,$Ns,$basisarsort,$Nstate,$tagssort)
@btime Hammap($instate,$zeta,$U,$Ns,$basisarsort,$Nstate,$tagssort)

```

---

<div class="post-metadata">

**Author:** ![eveningsilverfox](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eveningsilverfox/32/51932_2.png) [@eveningsilverfox](https://discourse.julialang.org/u/eveningsilverfox)\
**Post date:** [March 14, 2025, 10:21pm UTC](https://discourse.julialang.org/t/matrix-linear-map-function-exponential-multiplied-with-vector-in-exponentialutilities-jl/126957/10 "2025-03-14T22:21:03Z")

</div>

In this toy version, where I use the simple functions for creating the matrices, the map-based way is actually faster.

```julia
using MKL
using LinearAlgebra
using ArnoldiMethod
using SparseArrays
using MKLSparse
using LinearMaps    
using BenchmarkTools

function matmap(x,aa,bb)
    Nx = size(x,1);
    y = zeros(eltype(x), Nx);
    for gh = 1:Nx
        @views y[gh] += aa*x[gh];
        if gh>1 && gh<Nx
            @views y[gh] += -bb*x[gh-1]-bb*x[gh+1];
        end
        @views y[1] += -bb*x[2];
        @views y[Nx] += -bb*x[Nx-1];
    end

    return y
end

function matmat(aa,bb,N)
    rows = Float64[];
    cols = Float64[];
    vals = Float64[];
    for gh = 1:N
        # Mmat[gh,gh] = aa;
        append!(rows, gh);
        append!(cols, gh);
        append!(vals, aa);
    
        if gh > 1 && gh < N
            # Mmat[gh,gh-1] = Mmat[gh,gh+1] = -bb;
            append!(rows, gh);
            append!(cols, gh-1);
            append!(vals, -bb);
            append!(rows, gh);
            append!(cols, gh+1);
            append!(vals, -bb);
        end 
    end
    # Mmat[1,2] = Mmat[N,N-1] = -bb;
    append!(rows, 1);
    append!(cols, 2);
    append!(vals, -bb);
    append!(rows, N);
    append!(cols, N-1);
    append!(vals, -bb);

    # Mmatsparse = sparse(Mmat);
    Mmatsparse = sparse(rows,cols,vals,N,N);

    return Mmatsparse
end

function arnoldimap(N)
    decomp, history = partialschur(LinearMap(x->matmap(x,aa,bb), N), nev=1, which=:SR); # only solve for the ground state
    eval, evec = partialeigen(decomp);

    return eval, evec
end

function arnoldimat(Mat, N)
    decomp, history = partialschur(Mat, nev=1, which=:SR); # only solve for the ground state
    eval, evec = partialeigen(decomp);

    return eval, evec
end

aa = 0.9; bb = 0.5; N = 6000;
mat = matmat(aa,bb,N);

@btime arnoldimap($N);
@btime arnoldimat($mat, $N);

```

---

<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 14, 2025, 11:19pm UTC](https://discourse.julialang.org/t/matrix-linear-map-function-exponential-multiplied-with-vector-in-exponentialutilities-jl/126957/11 "2025-03-14T23:19:20Z")

</div>

> [@eveningsilverfox](#):
>
> Yes, I should use findmin instead. The partialsortperm is a vestigial code from many other stuff I have been copying from.

If you are doing linear search over and over you should probably look for a different data structure. Also, you are allocating a temporary array just for the argument to partialsortperm, which you can easily eliminate. Saving the index array alleviates both these problems, but might be suboptimal.

---

<div class="post-metadata">

**Author:** ![eveningsilverfox](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eveningsilverfox/32/51932_2.png) [@eveningsilverfox](https://discourse.julialang.org/u/eveningsilverfox)\
**Post date:** [March 15, 2025, 12:04am UTC](https://discourse.julialang.org/t/matrix-linear-map-function-exponential-multiplied-with-vector-in-exponentialutilities-jl/126957/12 "2025-03-15T00:04:07Z")

</div>

Yea, the map-based Arnoldi makes multiple (~100 it seems) calls to the massive array of basis elements (basisarsort) and the hash lookups with argmin (tagssort) every time it calls the function which creates the map (Hammap).  
This problem will persist even on using expv! from ExponentialUtilities.jl. Instead I am using a RK4 evolution, which I guess only preserves the norm of the vector upto some time.

While the matrix is sparse, it depends intricately on the ``content" of each basis element, and not on the index/tag/hash-output. So the hash lookup seems inevitable one way or the other.