Julia中自定义AbstractArray子类型矩阵运算精度偏差的原因及简化解决方法咨询

你遇到的精度问题根源很明确：你的自定义adjoint方法返回了FooArray类型，导致LinearAlgebra调用了通用的矩阵乘法实现，而原生矩阵R的转置乘法使用了经过BLAS优化的专用实现。这两种实现的浮点计算顺序、舍入策略不同，所以产生了微小的数值差异，累积后会变得明显。

不用手动重定义大量函数，有几种更简便的解决方案：

方案1：让adjoint返回原生矩阵的转置类型 #

既然你在转置时用了zeros作为新的vec（说明转置后的vec对你来说不重要），直接让adjoint返回底层数据的原生转置类型即可，这样运算会自动复用Julia优化后的乘法路径：

Base.adjoint(a::FooArray) = adjoint(a.data)

现在运行S'y - R'y就会得到全0向量，因为S'和R'现在都是原生的Adjoint矩阵类型，调用的是完全相同的BLAS乘法实现。

方案2：让LinearAlgebra识别你的底层数据 #

如果你需要保留FooArray的结构（比如后续还要访问列名称），可以定义Base.parent方法，让LinearAlgebra框架直接访问你的data字段，自动复用原生矩阵的优化操作：

Base.parent(fooarr::FooArray) = fooarr.data

这个方法会让很多LinearAlgebra函数（包括乘法、分解等）自动跳过自定义类型的包装，直接对底层数据进行运算，既保留了FooArray的元数据，又能获得原生的精度和性能。

额外建议：选择性保留自定义类型结构 #

如果某些运算需要返回FooArray（比如两个FooArray相乘后要保留新的列名称），你只需要针对性地重定义这些特定运算，而不是所有函数。例如：

function *(a::FooArray{T,2}, b::FooArray{S,2}) where {T,S}
    new_data = a.data * b.data
    # 这里可以根据需求生成新的vec，比如合并列名称或者自定义逻辑
    new_vec = # 你的列名称生成逻辑
    FooArray(new_data, new_vec)
end

对于像FooArray * AbstractVector这种返回普通向量的操作，完全不需要重定义，LinearAlgebra会自动使用底层数据的优化实现。

这样既解决了精度问题，又避免了维护大量函数的麻烦。

内容的提问来源于stack exchange，提问作者Claudio