实现矩阵函数的最有效方法？答案

【问题标题】：Most efficient way to implement matrix functions?实现矩阵函数的最有效方法？
【发布时间】：2017-12-06 03:41:53
【问题描述】：

我一直在尝试找出处理包含变量的二维数组（即矩阵）的最有效方法。我通过连接较小的数组来构造一个大矩阵，并用一个变量加权。作为一个简化的例子，这就是我目前所做的。

function myMatrix(x::Float64)
    M = vcat(hcat(x*A, C), hcat(C, 2*x*B))
    return M
end

const A = ones(2,2)  
const B = ones(2,2)
const C = zeros(2,2)
y_values = collect(0:0.01:10)

for y in y_values
    eivals, eivecs = eig(myMatrix(y))
    dosomething(eivals,eivecs) 
end

问题：我有更多（和更复杂的非对角线）矩阵，最终大小非常大。此刻，我调用了这个函数几百次。有没有人建议如何使这个过程在运行时更有效？

提前致谢

【问题讨论】：

撇开核心问题不谈，这将非常缓慢，因为您在 myMatrix 中使用了非常量全局变量，这是（我认为）性能提示中的第一个禁忌官方文档的部分。
谢谢，我改了。
我应该提一下，在我的原始程序中，矩阵实际上取决于代码的一些早期部分，特别是一些输入参数（每次程序运行都设置一次）。但是，我遵循了那里的性能提示并将它们声明为全局。
制作一个预分配版本，您可以在其中更改传入的某些输出变量的值，而不是每次都创建一个新矩阵。
我不确定我理解你的意思，你能具体说明一下吗？

标签： arrays performance julia

【解决方案1】：

编辑：根据下面的评论，这不能回答问题。

你的意思是构建矩阵吗？一种符合意识形态的方法是使用块矩阵：

M = [x*A   C
    C    2x*B]

不确定这在运行时是否最有效，抱歉。

【讨论】：

不，我可以构建矩阵（至少我认为）。我的问题是我做了几百次连接过程（即函数调用），这似乎有点白痴。我正在构建的矩阵是由几百个较小的矩阵构成的，而不仅仅是像这里的四个，所以我不能像你那样明确地写出来。我希望可能有一种计算效率更高的方法来处理这个矩阵作为 x 的函数，而不是基本上为每个 x 值重复构造矩阵。

【解决方案2】：

该解决方案使用 cmets 中某处提出的预分配，并不能直接解决我的问题（因为我构造矩阵的方式），但对于阅读本文的某些人可能仍然有用。此外，我不保证最初使用的“最有效”声明，因为这似乎在很大程度上取决于您的目的、矩阵的大小等。该方法在 Julia 的性能提示部分也有提及。

由于出现了一些混淆，请考虑以下示例：

function myMatrix(x::Float64)
    M = vcat(hcat(x*A, C), hcat(D, 2*x*B))
    return M
end

function doSomething(A::Array{Float64,2})
    nothing
end

const ArraySize = 1000
const A = ones(ArraySize,ArraySize)  
const B = ones(ArraySize,ArraySize)
const C = rand(ArraySize,ArraySize)
const D = rand(ArraySize,ArraySize)


for i = 1:1000
    ret = myMatrix( convert(Float64,i) )
    doSomething(ret)
end

这实际上什么都不做，只是根据一个参数从初始矩阵构造一个 BlockMatrix(function)。我在想这种重复构造是多余的，实际上，可以通过编写来为矩阵预分配内存：

function xinc!(ret::Array{T,2}, x::T) where T
    ret[1:ArraySize, 1:ArraySize] = x*A
    ret[1:ArraySize, ArraySize+1:2*ArraySize] = C
    ret[ArraySize+1:2*ArraySize, 1:ArraySize] = D
    ret[ArraySize+1:2*ArraySize, ArraySize+1:2*ArraySize] = 2*x*B
    nothing
end

function doSomething(A::Array{Float64,2})
    nothing
end

const ArraySize = 1000
const A = ones(ArraySize,ArraySize)  
const B = ones(ArraySize,ArraySize)
const C = rand(ArraySize,ArraySize)
const D = rand(ArraySize,ArraySize)



ret = Array{Float64}(2*ArraySize, 2*ArraySize)
for i = 1:1000
    xinc!(ret, convert(Float64,i))
    doSomething(ret)
end

对我来说，第二个代码在 9.866 秒内执行，而第一个代码占用 38.076 秒。

编辑：回应之前的评论，如果我写了

function xinc!(ret::Array{T,2}, x::T) where T
    ret = [x*A  C
            D   2*x*B]
    nothing
end

代码执行需要 16.173 秒。我不知道为什么，但是这种分配矩阵的方式要慢得多。

【讨论】：