expv slow for Bidiagonal and Tridiagonal

[sethaxen]

expv is very slow for Bidiagonal and Tridiagonal compared to their equivalent SparseMatrixCSC representations. A minimal example:

using LinearAlgebra, BenchmarkTools, SparseArrays
n = 1000
Atri = Tridiagonal(randn(n, n))
Atri_sp = sparse(Atri)
Abi = Bidiagonal(randn(n, n), :U)
Abi_sp = sparse(Abi)
b = randn(n)

julia> @btime expv(1.0, $Atri, $b);
  1.878 ms (101 allocations: 786.73 KiB)

julia> @btime expv(1.0, $Atri_sp, $b);
  318.413 μs (104 allocations: 849.47 KiB)

julia> @btime expv(1.0, $Abi, $b);
  1.889 ms (118 allocations: 929.56 KiB)

julia> @btime expv(1.0, $Abi_sp, $b);
  257.115 μs (95 allocations: 762.41 KiB)

The culprit is parameters. e.g.

julia> @btime $(ExponentialAction.parameters)(1.0, $Atri, 1, 50)
  1.684 ms (1 allocation: 896 bytes)
(m = 44, s = 1)

julia> @btime $(ExponentialAction.parameters)(1.0, $Atri_sp, 1, 50)
  4.893 μs (1 allocation: 896 bytes)

This partially is because Atri * Atri is a SparseMatrixCSC, which is then multiplied by Atri:

julia> @btime $(Atri^2) * $Atri;
  28.066 ms (27 allocations: 265.16 KiB)

julia> @btime $(Atri^2) * $Atri_sp;
  60.443 μs (6 allocations: 280.34 KiB)

Because there's no specialized overload for *(::SparseMatrixCSC, ::Tridiagonal), the product is very slow.

Also, opnorm(::SparseMatrixCSC, 1) has a fast custom overload, while opnorm(::Tridiagonal, 1) does not:

julia> @btime $(ExponentialAction.opnormest1)($Atri)
  1.684 ms (0 allocations: 0 bytes)
6.882369544243176

julia> @btime $(ExponentialAction.opnormest1)($Atri_sp)
  4.485 μs (0 allocations: 0 bytes)
6.882369544243176

A simple workaround for these matrices is a custom overload for parameters to sparsify the matrix before passing it to the standard parameters, which is cheap:

julia> @btime sparse($Atri);
  5.294 μs (5 allocations: 55.03 KiB)

This gets us down to

julia> @btime expv(1.0, $Atri, $b);
  199.531 μs (106 allocations: 841.77 KiB)

julia> @btime expv(1.0, $Abi, $b);
  202.301 μs (121 allocations: 969.06 KiB)