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

General Usage
, , ,
1

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

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

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.

2

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).

4 Likes
3

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?

4

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!.)

6

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.

7

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

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

8

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?)

9

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.

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

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

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);
11

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.

12

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.