在茱莉亚合并repmat和转置

Question

我将一些代码从R移植到julia以熟悉该语言，并且我发现了一些不能平滑转换的模式。考虑下面的函数，

# Ricatti-Bessel and derivatives up to nmax, vectorised over x 
function rb(x, nmax) 

    n = 1:nmax 
    nu = 0.5 + [0, n] 

    bj = hcat([besselj(nu, _x) for _x in x]...).' 
    #^first question^
    sq = repmat(sqrt(pi/2*x), 1, nmax+1) 

    bj .*= sq 
    xm = repmat(x, 1, nmax) 
    nm = repmat(n', length(x), 1) 
    #^second question^

    dpsi = bj[:,n] - nm .* bj[:,n+1] ./ xm 
    psi = bj[:,n+1] 
    return psi, dpsi # it'd be nice to return a "named list" instead 

end 
# rb(1:5,3) 

• 第一个问题：这是获得使用besselj() n最大列长度（X）列的矩阵的最佳方式？在找到一种可行的模式之前，我不得不挠挠头。

• 第二个问题：我发现自己经常不得不转换对象，在repmat中和/或在repmat之外，有没有其他方法可以指定输出大小和填充方向（按行还是在col-wise）？

也许我采取了错误的做法与整个事情：我使用矢量化功能（以R，和Matlab的老人回忆）的工作，因为他们通常是快速的最短路线线性代数的例程。 保持x是一个标量，并且只在最高级别循环会更有意义吗？我担心通过这样做，我将无法利用BLAS等的快速矩阵/矢量函数，并基本上将它们重写在茱莉亚中，更不用说可读性明显下降了。 我应该强调我对最佳性能感兴趣，因为这个函数将在很多次内部被调用，对于x的许多值。

Answer 1

关于第一个问题，我会跟下面的矩阵理解替换：

nu = (0:nmax)+0.5 
bj = [besselj(i,j) for j in x, i in nu] 

对于你的第二个，我认为一个好的原则，在朱莉娅编写高性能的代码是避免不必要的分配（和读取当然！）当Julia生成非常快速的指令时 - 这就是为什么for循环完美无缺，除矢量化线性代数（例如矩阵乘法）之外并不重要。它不擅长的工作是避免分配不需要的内存（像sq这样的临时矩阵）。我更换了bj/sq线，只是下面的

nu = (0:nmax)+0.5 
bj = [besselj(i,j)*sqrt(pi/2*j) for j in x, i in nu] 

这是很好的，因为它只是一个分配，在此之前我们已经（假设我们从我的回答Q1开始）：

1. 分配bj
2. 分配sq
3. 分配bj.*sq和重新绑定bj这种新的内存

（注意.*=是不就地操作！）

您的“命名名单”请求可能最好通过使type这个函数的返回（现在无论是满足这根本不是一个昂贵的操作，并且在例如非常普遍Julia的矩阵分解代码，其中需要返回多个值）。或者，您可以返回一个Dict，但这并不感觉习惯。

对于dpsi这一行，我会给你两个选项。首先是另一种矩阵理解：

dpsi = [ bj[i,j] - j * bj[i,j+1]/i 
      for i in 1:length(x), j in 1:nmax] 

，另一个是环Y：

dpsi = zeros(length(x),nmax) 
for i in 1:length(x), j in 1:nmax 
    dpsi[i,j] = bj[i,j] - j * bj[i,j+1]/i 
end 

在这两种情况下，我回避的临时的分配。同样，你的原单有如下分配：

1. 对于​​
2. 对于nm
3. 对于bj[:,n]（这将改变是在0.4的视图）
4. 对于bj[:,n+1]（同上）
5. 对于
6. 对于nm .* bj[:,n+1] ./ xm
7. 对于th Ë整个结果

两者我提出的版本只有一个配置，也很可能是更接近问题的原来的数学声明

我的最终版本是

function myrb(x, nmax) 
    bj = [ besselj(i,j)*sqrt(pi/2*j) 
       for j in x, i in (0:nmax)+0.5] 
    dpsi = [ bj[i,j] - j * bj[i,j+1]/i 
       for i in 1:length(x), j in 1:nmax] 
    psi = bj[:,2:nmax+1] 
    return psi, dpsi 
end 

我不对于besselj知之甚少，但我猜这是整个事件中最慢的部分，所以在这个特殊情况下，所有这些在速度上可能都没有多大关系。标杆这条微的情况下建议尽可能多：

# original 
elapsed time: 9.7578e-5 seconds (7176 bytes allocated) 
elapsed time: 7.2644e-5 seconds (7176 bytes allocated) 
elapsed time: 7.5709e-5 seconds (7176 bytes allocated) 
# revised 
elapsed time: 2.7536e-5 seconds (728 bytes allocated) 
elapsed time: 2.7097e-5 seconds (728 bytes allocated) 
elapsed time: 1.6601e-5 seconds (728 bytes allocated) 

您可以用探查证实了这一点（虽然我不得不用更大的投入:)

@profile myrb(1:500,300) 
Profile.print() 

在我的机器有在采集样本429功能，其中426个位于Julia内部的bessel.jl文件中，2个位于dpsi，1个位于psi。