为什么 AbstractArray 的子类型会导致 Julia 中的矩阵运算不精确？

Question

我目前正致力于在 Julia 中创建 AbstractArray 的子类型，它允许您在数组本身之外存储向量。您可以将其视为“名称”列，元素类型作为 AbstractFloat 的子类型。因此，它与 NamedArray.jl 包有一些相似之处，但仅限于使用浮点数分配列（在矩阵的情况下）。

到目前为止我创建的结构（按照 guide 创建 AbstractArray 的子类型）定义如下：

struct FooArray{T, N, AT, VT} <: AbstractArray{T, N}
    data::AT
    vec::VT
    function FooArray(data::AbstractArray{T1, N}, vec::AbstractVector{T2}) where {T1 <: AbstractFloat, T2 <: AbstractFloat, N}
        length(vec) == size(data, 2) || error("Inconsistent dimensions")
        new{T1, N, typeof(data), typeof(vec)}(data, vec)
    end
end
@inline Base.@propagate_inbounds Base.getindex(fooarr::FooArray, i::Int) = getindex(fooarr.data, i)
@inline Base.@propagate_inbounds Base.getindex(fooarr::FooArray, I::Vararg{Int, 2}) = getindex(fooarr.data, I...)
@inline Base.@propagate_inbounds Base.size(fooarr::FooArray) = size(fooarr.data)
Base.IndexStyle(::Type{<:FooArray}) = IndexLinear()

这似乎足以创建 fooArray 类型的对象并用它做一些简单的数学运算。但是，我观察到一些基本函数（例如矩阵向量乘法）似乎不精确。例如，下面应该始终 return 一个 0.0 的向量，但是：

R = rand(100, 3)
S = FooArray(R, collect(1.0:3.0))
y = rand(100)
S'y - R'y
3-element Vector{Float64}:
 -7.105427357601002e-15
 0.0
 3.552713678800501e-15

虽然差异非常小，但它们可以通过许多不同的计算快速累加，从而导致重大错误。

这些差异从何而来？

通过宏 @code_llvm 查看计算表明使用了与 LinearAlgebra 明显不同的 matmul 函数（有其他细微差别）：

@code_llvm S'y
   ...
   @ C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\matmul.jl:111 within `*'
   ...
@code_llvm S'y
   ...
   @ C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\LinearAlgebra\src\matmul.jl:106 within `*'
   ...

在我们的 FooArray 对象上重新定义 adjoint 和 * 函数提供了预期的正确结果：

import Base: *, adjoint, /
Base.adjoint(a::FooArray) = FooArray(a.data', zeros(size(a.data, 1)))
*(a::FooArray{T, 2, AT, VT} where {AT, VT}, b::AbstractVector{S}) where {T, S} = a.data * b
S'y - R'y
3-element Vector{Float64}:
 0.0
 0.0
 0.0

但是，此解决方案（也在 NamedArrays here 中完成）需要定义和维护各种函数，而不仅仅是 base 中的标准函数，添加更多和更多的依赖只是因为这个小误差范围。

有没有更简单的方法来解决这个问题，而无需重新定义每个操作以及其他包中可能的许多其他函数？

我在 Windows 64 位系统上使用 Julia 1.6.1 版。

Answer 1

是的，矩阵乘法的实现会因您的数组类型而异。内置 Array 将使用 BLAS，而您的自定义 fooArray 将使用通用实现，并且由于浮点运算的非关联性，这些不同的方法确实会产生不同的值——请注意它们可能与基本事实不同，即使对于内置 Arrays!

julia> using Random; Random.seed!(0); R = rand(100, 3); y = rand(100);

julia> R'y - Float64.(big.(R)'big.(y))
3-element Vector{Float64}:
 -3.552713678800501e-15
  0.0
  0.0

您可以将自定义数组实现为 DenseArray，这将确保它使用相同的（支持 BLAS 的）代码路径。你只是 need to implement a few more methods，最重要的是 strides 和 unsafe_convert:

julia> struct FooArray{T, N} <: DenseArray{T, N}
          data::Array{T, N}
       end
       Base.getindex(fooarr::FooArray, i::Int) = fooarr.data[i]
       Base.size(fooarr::FooArray) = size(fooarr.data)
       Base.IndexStyle(::Type{<:FooArray}) = IndexLinear()
       Base.strides(fooarr::FooArray) = strides(fooarr.data)
       Base.unsafe_convert(P::Type{Ptr{T}}, fooarr::FooArray{T}) where {T} = Base.unsafe_convert(P, fooarr.data)

julia> R = rand(100, 3); S = FooArray(R); y = rand(100)
       R'y - S'y
3-element Vector{Float64}:
 0.0
 0.0
 0.0

julia> R = rand(100, 1000); S = FooArray(R); y = rand(100)
       R'y == S'y
true

为什么 AbstractArray 的子类型会导致 Julia 中的矩阵运算不精确？

Why does a subtype of AbstractArray result in imprecise matrix operations in Julia?

arrays

floating-accuracy

julia