Large Matrix Takes Forever

General Usage Performance
1

Hi, I am new to Julia and I am writing a program that makes a random matrix drawn from the Haar distribution. It works fine for small matrices ~200x200 but when I try for 1000x1000 the program takes a very long time to run. I am not sure what the issue is here.

using Random, LinearAlgebra
using Plots

function mkHaar(dimMat)
    #initialize random number generatore and create complex 4x4 normal random matrix
    rng = MersenneTwister()
    rndMat = randn(rng, ComplexF64, (dimMat, dimMat))

    Z = qr(rndMat)
    L = Diagonal(Z.R)/abs.(Diagonal(Z.R))

    rndHaar = Z.Q*L

    #We can check to see that rndHaar is Unitary
    #println(adjoint(rndHaar)*rndHaar)

    #and that the eigenvalue phases are uniformly distributed
    EV = eigvals(rndHaar)
    phaseEV = angle.(EV)

    #Make histogram
    gr(size = (300,300), legend = false)
    #create binning breaks
    binbrks = range(-pi, stop = pi, length = 25)
    p = histogram(phaseEV,bins=binbrks)
    display(p)

end
2

It looks like matrix multiplication with Z.Q (which is not actually a plain matrix but a special type for the result of the QR factorization) is very slow because it ends up calling the getindex function here: julia/qr.jl at 8639dc1dbfb9d0a8df50e1eb95e68cc6be880bd3 Ā· JuliaLang/julia Ā· GitHub which looks pretty slow to me.

You can get around this by just copying it to a matrix first:

rndHaar = Matrix(Z.Q) * L

Thereā€™s probably a better way to fix this in Julia itself, but I donā€™t know that part of the code well enough.

3 Likes
3

Hereā€™s an implementation that will work (you can remove the two Compat. prefixes if you are on 1.0):

4

There is a fast method for multiplying the special ā€œQā€ matrix type by a StridedMatrix: https://github.com/JuliaLang/julia/blob/master/stdlib/LinearAlgebra/src/qr.jl#L567
but this doesnā€™t get called for multiplying by a Diagonal matrix, for which it ends up calling *(::AbstractMatrix, ::Diagonal) instead.

To avoid having an explosion of the number of *(...) methods we need, it might be nice to have some kind of ExpensiveIndexing trait that generic methods can check to see if they need to make a Matrix copy before executing indexing-heavy algorithms.

3 Likes
5

Hi,

Thanks for all the replies they were very helpful. For now I think I will just add the Matrix() since it is the simplest but Iā€™ll do some performance checks and see whatā€™s the best way to handle this in the future.