Yeppp-package

New to Julia
,
1

Hi, I have seen the package Yeppp been mentioned here (Fast logsumexp). This might be an ignorant question but when adding Yeppp it wont install, and I get some errors telling me to Pkg.build instead, which also does not work, and I wonder why.
Is it no longer preferred to use the package? Or maybe the performance is now included in LinearAlgebra?

New to Julia, using 1.8.3, and just want to check if using Yeppp makes a difference in improving parts of a numpy/scipy-based-computation , specifically the step sum(log.(eigvals(Symmetric(A)))) where A is Matrix of Float64, different sizes every computation ranging from 2x2 to 400x400

Best,

2

The repo for that package has been archived and has had no commits for four years, which is quite long in Julia time. I don’t know if it’s been superceded by something or just stopped being developed.

By the way, a simple improvement is to use sum(log, x) rather than sum(log.(x)). It should make a tiny difference compared to all those eigs, but a difference nonetheless.

2 Likes
4

I’d recommend LoopVectorization.jl for vectorizing functions like log.
But I’d be surprised if even base log takes much time compared to eigvals.

1 Like
5

@gustaphe @Elrod Thank you for the suggestions.

I see sum(log, ex) is better, at least it allocates less.
Can I ask a follow-up question regarding eigvals? I wonder if there is an obvious code that could improve performance that I am missing. Would LoopVectorization.jl be something to look into for this operation too?

6

If you want better performance of the log sum:

julia> using LoopVectorization

julia> x = rand(511);

julia> @btime vmapreduce(log, +, $x)
  277.449 ns (0 allocations: 0 bytes)
-487.3246594217604

julia> @btime sum(log, $x)
  2.490 μs (0 allocations: 0 bytes)
-487.3246594217602

Unfortunately, it won’t help eigvals just yet.

5 Likes
7

Wow, that’s almost 10x!

1 Like
8

If you are dealing with very large sparse matrices A you might be interested in O(N) spectral methods, like the Kernel Polynomial Method and relatives. See e.g. Rev. Mod. Phys. 78, 275 (2006) - The kernel polynomial method

The basic idea is to express the function of A that you want to trace over (in this case log) as an expansion in Chebyshev polynomials P_n(x), which end up giving you the sum you want in terms of psi'*P_n(A)*psi, where psi is a random vector. These methods are crazy efficient for very large sparse matrices.

Probably for 400x400 dense matrices this doesn’t make much sense though.

2 Likes
9

@pablosanjose @Elrod Thank you, I just tried vmapreduce() and it helps in this calculation, it seems awesome. As you mentioned earlier, it makes a tiny difference compared to the eigvals!(Hermitian(A)). Do you know if it is possible to improve this eigvals-computation?
I forgot to mention the matrices are dense, otherwise I would try out the O(N) spectral method you mentioned above. Thanks alot for the suggestion.

10

If you have an Intel or even an AMD CPU, using MKL should help.

1 Like
11

Note that this is the same thing as the log-determinant logdet(A), which is much faster than computing eigenvalues. Since it seems like you know that your matrix A is positive definite (as otherwise log(λ) will throw an error), you can often do even better still by using logdet(cholesky(Symmetric(A))):

julia> using LinearAlgebra, BenchmarkTools

julia> BLAS.set_num_threads(1) # just a single core for benchmarking

julia> A = randn(1000,1000); A = A'A; # random SPD matrix;

julia> @btime sum(log, eigvals(Symmetric($A)))
  68.893 ms (11 allocations: 7.99 MiB)
5902.6715053545195

julia> @btime logdet($A)
  18.700 ms (3 allocations: 7.64 MiB)
5902.671505354289

julia> @btime logdet(cholesky(Symmetric($A)))
  11.704 ms (2 allocations: 7.63 MiB)
5902.671505354335

The log-determinant is a useful function with lots of beautiful properties!

12 Likes