在 julia 中用 3 个向量和 n 个向量广播向量乘法

Question

我有一个 3 向量 c = [0.7, 0.5, 0.2]，我想将它与 n 向量 x = rand((-1,1),n) 中的所有内容相乘，这样我得到一个结果 n+2 向量 y，其中 y[i] == x[i]*c[3] + x[i-1]*c[2] + x[i-2]*c[1]

我应该如何在 julia 中执行此操作？我觉得应该有一种方法可以将较小的 3 向量广播到 n 向量中的所有值。对于边缘情况，如果 i-1 或 i-2 超出范围，我只希望这些组件为零。

Answer 1

这是一种使用 ShiftedArrays.jl.

的方法

using ShiftedArrays

c = [0.7, 0.5, 0.2]

创建 x 的滞后版本，初始为零：

x = 1:5
xminus1 = lag(x, 1, default=0)
xminus2 = lag(x, 2, default=0)

水平连接向量并使用矩阵乘法 c:

X = [xminus2 xminus1 x]
X * c

这是 X 和 X * c 在 REPL 中的样子：

julia> X = [xminus2 xminus1 x]
5×3 Array{Int64,2}:
 0  0  1
 0  1  2
 1  2  3
 2  3  4
 3  4  5

julia> X * c
5-element Array{Float64,1}:
 0.2
 0.9
 2.3
 3.7
 5.1

请注意，这会生成长度为 length(x) 的输出向量，而不是 length(x) + 2。我不确定输出的长度如何有意义 length(x) + 2，正如您在问题中所要求的那样。

Answer 2

如果我正确理解你的问题，你需要一个卷积，在标准卷积中，向量 c 会被反转。你可以使用例如DSP.jl为此。

这是你想要的吗？

julia> using DSP

julia> c = [0.7, 0.5, 0.2]
3-element Array{Float64,1}:
 0.7
 0.5
 0.2

julia> conv([10, 100, 1000, 10000], reverse(c))
6-element Array{Float64,1}:
    1.9999999999996967
   25.0
  257.0000000000003
 2569.9999999999995
 5700.0
 6999.999999999998

您也可以使用 LinearAlgebra 模块中的 dot 手动实现它，如下所示：

julia> using LinearAlgebra

julia> x = [10, 100, 1000, 10000]
4-element Array{Int64,1}:
    10
   100
  1000
 10000

julia> y = [0;0;x;0;0]
8-element Array{Int64,1}:
     0
     0
    10
   100
  1000
 10000
     0
     0

julia> [dot(@view(y[i:i+2]), c) for i in 1:length(x)+2]
6-element Array{Float64,1}:
    2.0
   25.0
  257.0
 2570.0
 5700.0
 7000.0

Answer 3

我有一个包可以做这些事情。最简单的用法是这样的：

julia> c = [0.7, 0.5, 0.2]; # from question

julia> x = [10, 100, 1000, 10_000]; # from another answer

julia> using Tullio, OffsetArrays

julia> @tullio y[i] := x[i]*c[3] + x[i-1]*c[2] + x[i-2]*c[1]
2-element OffsetArray(::Vector{Float64}, 3:4) with eltype Float64 with indices 3:4:
  257.0
 2570.0

julia> @tullio y[i] := x[i+k-3] * c[k] # sum over all k, range of i that's safe
2-element OffsetArray(::Array{Float64,1}, 3:4) with eltype Float64 with indices 3:4:
  257.0
 2570.0

因为 eachindex(c) == 1:3，这是求和的 k 值的范围，i 的范围尽可能大，因此 i+k-3 留在 eachindex(x) == 1:4.

里面

要通过在每个方向用两个零填充 x 来扩展 i 的范围，请写入 pad(i+k-3, 2)。要计算生成普通的基于 1 的数组所需的 i 的移位，请在左侧写入 i+_（然后 -3 没有区别）。那么：

julia> @tullio y[i+_] := x[pad(i+k, 2)] * c[k]
6-element Array{Float64,1}:
    2.0
   25.0
  257.0
 2570.0
 5700.0
 7000.0

在较大的数组上，这不会很快（目前），因为它必须在每一步检查它是在 x 内还是在填充内。 DSP.conv 很可能在这方面更聪明一些。（编辑 - DSP.jl 似乎永远不会比这个 c 更快；对于长度为 1000 的内核，x 中有 10^6 个元素会更快。）

在 julia 中用 3 个向量和 n 个向量广播向量乘法

broadcasting vector multiplication with a 3 vector and n vector in julia

vector

linear-algebra

julia