在朱莉娅结合repmat和transpose

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)

• 第一个问题:这是获得使用nmax列和长度(x)行的矩阵的最佳方法besselj()吗？在找到一个有效的模式之前,我不得不抓了一会儿.

• 第二个问题:我发现自己经常不得不转换物品,在里面和/或外面重新定位,是否有另一种选择,我可以指定输出大小和填充方向(行方式或列方式)？

也许我对整个事情采取了错误的方法:我习惯于使用矢量化函数(在R中,以及Matlab的旧记忆中),因为它们通常是线性代数快速例程的最短路径.保持xa标量是否更有意义,并且仅在最高级别循环？我担心通过这样做,我将无法利用BLAS等的快速矩阵/向量函数,并且基本上在julia中重写它们,更不用说明显的可读性损失了.我应该强调我对最佳性能感兴趣,因为对于x的许多值,这个函数将被内部调用很多次.

Answer 1

对于你的第一个问题,我将用以下矩阵理解代替它:

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

对于你的第二个,我认为在Julia中编写高性能代码的一个好原则是避免不必要的分配(当然,阅读性能提示!)Julia生成非常快速的指令 - 这就是为什么for循环非常好,并且它除了线性代数(例如矩阵乘法)之外,对矢量化其他任何东西都不是很关键.它做得不好的是避免分配不必要的内存(像你这样的临时矩阵sq).我用以下代替了bj/ sqlines

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. 对于 xm
2. 对于 nm
3. 对于bj[:,n](这将改为0.4中的视图)
4. 对于bj[:,n+1](同上)
5. 对于 nm .* bj[:,n+1]
6. 对于 nm .* bj[:,n+1] ./ xm
7. 对于整个结果

我建议的两个版本只有一个分配,可能更接近问题的原始数学陈述

我的最终版本是

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个bessel.jl在Julia内部,2个用于dpsi,1个用于psi.