# Large Matrix Takes Forever

**URL:** https://discourse.julialang.org/t/large-matrix-takes-forever/21328
**Category:** Performance
**Tags:** linearalgebra
**Created:** [March 1, 2019, 3:23am UTC](https://discourse.julialang.org/t/large-matrix-takes-forever/21328 "2019-03-01T03:23:18Z")
**Posts on this page:** 5
**Page:** 1

<div class="post-metadata">

### Author: ![olabaz](https://avatars.discourse-cdn.com/v4/letter/o/d07c76/32.png) [@olabaz](https://discourse.julialang.org/u/olabaz)
#### Post date: [March 1, 2019, 3:23am UTC](https://discourse.julialang.org/t/large-matrix-takes-forever/21328/1 "2019-03-01T03:23:18Z")

</div>

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.

```julia
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

```

---

<div class="post-metadata">

### Author: ![rdeits](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rdeits/32/286_2.png) [@rdeits](https://discourse.julialang.org/u/rdeits)
#### Post date: [March 1, 2019, 4:17am UTC](https://discourse.julialang.org/t/large-matrix-takes-forever/21328/2 "2019-03-01T04:17:28Z")

</div>

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](https://github.com/JuliaLang/julia/blob/8639dc1dbfb9d0a8df50e1eb95e68cc6be880bd3/stdlib/LinearAlgebra/src/qr.jl#L515-L521) which looks pretty slow to me.

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

```julia
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.

---

<div class="post-metadata">

### Author: ![jebej](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jebej/32/1784_2.png) [@jebej](https://discourse.julialang.org/u/jebej)
#### Post date: [March 1, 2019, 10:56pm UTC](https://discourse.julialang.org/t/large-matrix-takes-forever/21328/3 "2019-03-01T22:56:21Z")

</div>

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

> <https://github.com/jebej/Schrodinger.jl/blob/d8ff0197b0f01bc4634edb7548e7c7a18194ffb7/src/library/random.jl#L18>

---

<div class="post-metadata">

### Author: ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [March 2, 2019, 4:34am UTC](https://discourse.julialang.org/t/large-matrix-takes-forever/21328/4 "2019-03-02T04:34:42Z")

</div>

> [@rdeits](#):
>
> There’s probably a better way to fix this in Julia itself, but I don’t know that part of the code well enough.

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](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.

---

<div class="post-metadata">

### Author: ![olabaz](https://avatars.discourse-cdn.com/v4/letter/o/d07c76/32.png) [@olabaz](https://discourse.julialang.org/u/olabaz)
#### Post date: [March 4, 2019, 5:06pm UTC](https://discourse.julialang.org/t/large-matrix-takes-forever/21328/5 "2019-03-04T17:06:01Z")

</div>

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.
