Checking for code inefficiencies · helpdesk (published)

Hey everyone, I'm still struggling on spotting code inefficiencies and how to improve them.
Take these two functions below. Is there some inefficiency in them that could somehow be improved?

function w2gaussiancost(μa,Σa,μb, Σb)
    w2   = norm(μa - μb)^2 + tr(Σa + Σb - 2*(Σa^(1/2)* Σb*Σa^(1/2))^(1/2))
    return w2
end

function w2gaussianmap(μa,Σa,μb, Σb)
    A    = Σa^(-1/2)*(Σa^(1/2)* Σb*Σa^(1/2))^(1/2)*Σa^(-1/2)
    T(x) = μb + (A - μa)
    return T
end

Brenhin Keller (May 15 2021 at 02:13):

Davi Sales Barreira (May 15 2021 at 02:16):

Brenhin Keller (May 15 2021 at 02:19):

Gotcha. Yeah, you're pretty much dependent on the efficiency of the underlying lin alg methods in that case.
I don't see anything super obvious, but perhaps someone else will.

Brenhin Keller (May 15 2021 at 02:20):

If you suspect something isn't as efficient as it should be, you could try profiling and see if anyhing looks odd there

Brenhin Keller (May 15 2021 at 02:23):

I guess one other thing could be checking whether or not ^(1/2) is going to be the most efficient way to take the matrix square root

Pablo Zubieta (May 15 2021 at 02:43):

You might be able to avoid taking the square root of Σa twice by assigning its value to another variable and passing that to the expression

Pablo Zubieta (May 15 2021 at 02:44):

Brenhin Keller (May 15 2021 at 02:56):

Oh true yeah, the compiler almost certainly isn't going to be smart enough to do that for you

Michael Abbott (May 15 2021 at 03:10):

julia> begin
           n = 10
           Σa, Σb = (rand(n,n) for _ in 1:2)
           μa, μb = (rand(n) for _ in 1:2)
           rtΣa = @btime sqrt($Σa)
           @btime $rtΣa * $Σb * $rtΣa
           @btime w2gaussiancost($μa,$Σa,$μb,$Σb)  # original
           @btime w2gaussiancost_2($μa,$Σa,$μb,$Σb) # calls sqrt twice
       end;
  114.958 μs (10 allocations: 16.33 KiB)
  1.271 μs (2 allocations: 3.53 KiB)
  334.083 μs (49 allocations: 66.06 KiB)
  215.042 μs (25 allocations: 36.42 KiB)

Mason Protter (May 15 2021 at 05:30):

Mason Protter (May 15 2021 at 05:31):

Davi Sales Barreira (May 27 2021 at 20:27):

Davi Sales Barreira (May 27 2021 at 20:28):

Mason Protter (May 27 2021 at 21:05):

Okay, are you sure you need the actual inverse square root? Would a Cholesky factor work for you instead?

Mason Protter (May 27 2021 at 21:06):

Mason Protter (May 27 2021 at 21:11):

If you can't use Cholesky for some reason, I'd use either an SVD or Eigen decomposition and then use that to accellerate products of the form A * Σ^(-1/2)

Davi Sales Barreira (May 27 2021 at 21:54):

I think I do actually need it. But I'll see if perhaps the Cholesky decomp. works. Thanks for the input.

Mason Protter (May 27 2021 at 22:28):

@Davi Sales Barreira are you sure that your matrix is positive semi-definite, and not fully positive definite? For a semi-definite matrix, this whole thing is just ill conditioned when you take an inverse square root.

Davi Sales Barreira (May 27 2021 at 22:39):

It's a covariance matrix for the Multivariate Gaussian, so I think it is indeed positive semi-definite, but it is also Symmetric, so this might solve the ill condition. Doesn't it?

Davi Sales Barreira (May 27 2021 at 22:40):

 A    = Σa^(-1/2)*(Σa^(1/2)* Σb*Σa^(1/2))^(1/2)*Σa^(-1/2)

Which is the solution for the Optimal Transport map between multivariate gaussians. So I do think it should not be ill conditioned in any way, because this map is proved to always exist.

Mason Protter (May 27 2021 at 22:48):

Σa^(-1/2) is ill defined if the matrix is semi-definite, regardless of whether or not it's symmetric.

Mason Protter (May 27 2021 at 22:58):

Hmm, maybe it's okay actually, the square root does seem to improve it's conditioning quite a bit.

Davi Sales Barreira (May 27 2021 at 23:07):

I cannot prove it, but this result is very well established. So it's probably possible to prove it's indeed well conditioned.

Stream: helpdesk (published)

Topic: Checking for code inefficiencies

Davi Sales Barreira (May 15 2021 at 01:55):