超30k阶稠密非方阵并行最小二乘求解方案咨询

现有稠密非结构化矩阵满足等式 RG = A，各矩阵维度与数据类型如下：

• R：30k×40k 的 Float32 类型矩阵
• G：40k×50k 的 0/1 值矩阵（0、1元素占比大致相当）
• A：30k×50k 的 Float32 类型矩阵

已知 A 与 G，需要求解 R 的最小二乘解。可用计算资源包括数百个CPU核心、数百GB内存以及1张A40 GPU，示例代码基于Julia 1.7编写，也可接受其他可行技术方案。
第一个核心疑问：G 矩阵元素仅为0和1的特性是否可以被利用以优化求解效率？

基于Julia LinearAlgebra的多CPU方案测试 #

共测试了两种求解方法：Penrose inverse法与矩阵右除法，测试代码如下：

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')'

测试过程中发现两个问题：

• 线程数上限问题：执行 BLAS.set_num_threads(64) 或设置更高线程数时，BLAS.get_num_threads() 的返回值最高仅为32，存在线程数上限；

32个BLAS线程的运算速度反而慢于8线程 #

例如对维度为(4000, 9800) / (8500, 9800)的矩阵执行右除运算，8线程耗时不足50秒，32线程耗时超过55秒，经多次测试已排除编译耗时影响，不清楚该现象的成因，也不知道如何充分利用手上近100个CPU核心的算力。

测试发现矩阵右除法的运算速度快于Penrose inverse法，不清楚该现象是否符合预期：根据官方文档，左除运算符 \ 采用带主元选取的QR分解实现，但未查询到 pinv 与右除运算符 / 的具体实现算法（推测与 \ 算法一致，仅对矩阵做转置处理）。目前两种方法在目标规模的大矩阵上运算耗时极长，无法落地使用，希望了解如何调度多核CPU资源加速求解。

GPU方案测试情况 #

尝试使用 CUDA.jl 调用GPU加速求解，测试发现10k左右维度的矩阵可正常运行，pinv 运算耗时约1分钟，测试代码如下：

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)

但测试20k左右维度的矩阵时，程序直接抛出错误 InexactError: trunc(Int32, 4811456640)，推测是CUBLAS使用int32做索引导致（补充说明：报错原因是数组字节大小无法用31位整数表示）；对 CuArray 类型矩阵执行右除运算时，会抛出"DimensionMismatch("LU factored matrix A must be square!")的维度不匹配错误，不清楚需要手动选择哪种求解算法，也不确定对应算法在目标规模矩阵上是否仍会崩溃。

目前来看在Julia环境下无法便捷调用GPU解决该问题，希望得到指导：应该继续尝试适配GPU求解，还是采用多CPU方案更合适？

该问题为实际业务需求，请无需回复“没有人会需要这么大规模的最小二乘求解”类内容。

内容的提问来源于stack exchange，提问作者Adomas Baliuka