title: Parallellize least squares for large (> 30k x 30k) non-square dense matrices
body: <p>Let <code>RG = A</code> for dense unstructured matrices with shapes (e.g. roughly) <code>R</code>: (30k x 40k, entries float32) and <code>G</code>: (40k x 50k, entries either 0.0 or 1.0, roughly equally often) and of course <code>A</code>: (30k x 50k, entries float32).</p>
  <p>Given <code>A</code> and <code>G</code>, I want to find the least squares solution for <code>R</code>.</p>
  <p>I can use hundreds of CPU cores, hundreds of GB of RAM and also an A40 GPU. What is the best way to use such resources to solve the problem? <strong>I'm using Julia 1.7 in the examples below but I'm open to other options!</strong></p>
  <p><strong>First question: Can I somehow exploit that the entries of <code>G</code> are only zeros and ones?</strong></p>
  <h2>Trying to use Julia <code>LinearAlgebra</code> with many CPUs</h2>
  <p>I've tried two methods: &quot;Penrose inverse&quot; and &quot;right division&quot;</p>
  <pre><code>using LinearAlgebra
  @show BLAS.get_num_threads()
  # defaults to 8. Can change using BLAS.set_num_threads(N)
  
  # build toy problem (order of magnitude smaller sizes)
  R_true = rand(Float32, 3_000, 4_000) 
  G = rand([0., 1.], 4_000, 5_000)
  # note: using true/false here gives same results but is much slower!
  A = R_true * G
  
  # solve toy problem using matrix (right) division
  R_fitted_rdiv = A / G
  
  # solve toy problem using Penrose inverse
  R_fitted_pinv = (pinv(G') * A')'
  </code></pre>
  <p>First, setting <code>BLAS.set_num_threads(64)</code> (or any bigger number) actually only gives me <code>BLAS.get_num_threads()</code> returning 32. Apparantly that's an upper limit. Second,</p>
  <h3>using 32 BLAS threads is actually <em>slower</em> than using 8.</h3>
  <p>(e.g. performing right division with sizes (4000, 9800) / (8500, 9800) takes less than 50 seconds on 8 threads but more than 55 seconds on 32 threads. I ran things multiple times to exclude compilation time issues.) I don't know why this is or if it's normal. <strong>How can I make use of my computing power for this problem?</strong></p>
  <p><strong>I think that the matrix division is faster than the Penrose inverse method.</strong> Should this be expected? I don't know what either of the functions do exactly for these inputs. The docs say that left division (<code>\</code>) uses <a href="https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#Base.:-Tuple%7BAbstractMatrix,%20AbstractVecOrMat%7D" rel="nofollow noreferrer">pivoted QR factorization</a>. I couldn't find what algorithm(s) are used for <code>pinv</code> or right division (<code>/</code>) (although it's probably the same as <code>\</code> since they are related by transposing the matrices). I'd rather not delve too deeply because my knowledge in numerical linear algebra is quite limited.</p>
  <p><strong>The issue is that for my large matrices either method takes forever. Is there a way to make use of my ~100 cores somehow?</strong></p>
  <h2>Trying to use the GPU:</h2>
  <p>Using <code>CUDA.jl</code>, Matrices of size around 10k work fine and take a minute to <code>pinv</code>:</p>
  <pre><code>using CUDA
  @time matrix = CUDA.rand(Float32, 10_000, 10_500) # 0.003037 seconds (5 allocations: 160 bytes)
  @time pinv(matrix) #  57.417559 seconds (678 allocations: 172.094 KiB)
  </code></pre>
  <p>However, when I try to do matrices around size 20k, I get <strong>right away</strong> the error <code>InexactError: trunc(Int32, 4811456640)</code>. I assume this is due to <a href="https://github.com/JuliaGPU/CUDA.jl/issues/598" rel="nofollow noreferrer">CUBLAS using int32 for indexing</a>, even though I don't understand why it leads to an error in this case. <em>(edit: it's about the size of the array in bytes fitting into 31 bits.)</em></p>
  <p>Trying to use right division with <code>CuArray</code>s gives the error &quot;DimensionMismatch(&quot;LU factored matrix A must be square!&quot;)&quot;. I guess I have to choose a different algorithm manually? <strong>I don't know what it's called.</strong> (Although, it probably would still crash for large matrices...?)</p>
  <p>To summarize, it doesn't look like I can use the GPU from Julia easily to solve my problem. <strong>Should I keep trying to use the GPU for this task or stick to the many CPUs?</strong></p>
  <p><em>Yes this is really my problem, please refrain from commenting &quot;nobody should ever need such large least squares&quot;</em></p>