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

Question

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

到目前为止，我创建的结构（按照创建 的子类型的指南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并用它做一些简单的数学运算。然而，我发现一些基本函数（例如矩阵向量乘法）似乎不精确。例如，以下内容应始终返回 的向量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表明，使用了明显不同的matmul函数LinearAlgebra（还有其他细微差别）：

@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 `*'
   ...

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

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 此处完成）需要定义和维护各种函数，而不仅仅是 中的标准函数base，仅仅因为误差幅度很小，就添加了越来越多的依赖项。

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

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

Answer 1

是的，矩阵乘法的实现将根据您的数组类型而有所不同。内置函数Array将使用 BLAS，而您的fooArray自定义函数将使用通用实现，并且由于浮点运算的非关联性，这些不同的方法确实会产生不同的值 \xe2\x80\x94 并注意它们可能与基本事实，即使是内置的Array！

\n

julia> using Random; Random.seed!(0); R = rand(100, 3); y = rand(100);\n\njulia> R\'y - Float64.(big.(R)\'big.(y))\n3-element Vector{Float64}:\n -3.552713678800501e-15\n  0.0\n  0.0\n

\n

您可以将自定义数组实现为 a DenseArray，这将确保它使用相同的（启用 BLAS 的）代码路径。您只需要实现更多方法，最重要的strides是unsafe_convert：

\n

julia> struct FooArray{T, N} <: DenseArray{T, N}\n          data::Array{T, N}\n       end\n       Base.getindex(fooarr::FooArray, i::Int) = fooarr.data[i]\n       Base.size(fooarr::FooArray) = size(fooarr.data)\n       Base.IndexStyle(::Type{<:FooArray}) = IndexLinear()\n       Base.strides(fooarr::FooArray) = strides(fooarr.data)\n       Base.unsafe_convert(P::Type{Ptr{T}}, fooarr::FooArray{T}) where {T} = Base.unsafe_convert(P, fooarr.data)\n\njulia> R = rand(100, 3); S = FooArray(R); y = rand(100)\n       R\'y - S\'y\n3-element Vector{Float64}:\n 0.0\n 0.0\n 0.0\n\njulia> R = rand(100, 1000); S = FooArray(R); y = rand(100)\n       R\'y == S\'y\ntrue\n

\n