如何计算将矩阵保存为向量的矩阵乘法答案

【问题标题】：How to calculate matrix multiplication in which matrix is saved as vector如何计算将矩阵保存为向量的矩阵乘法
【发布时间】：2020-09-21 02:08:53
【问题描述】：

我有两个对称矩阵A、B 和一个向量X。 A的维数是n×n，B的维数是n×n，X的维数是n×1。让矩阵的ith行和jth列的元素A 表示为A[i,j]。

由于A是对称的，所以只保存A的上三角矩阵的每一列。矩阵A保存为数组：

Vector_A = [A[1,1],
            A[1,2], A[2,2],
            A[1,3], A[2,3], A[3,3],
            A[1,4], A[2,4], A[3,4], A[4,4],
            ...,
            A[1,n], A[2,n], ..., A[n,n]]

矩阵B 以与矩阵A 相同的格式保存。现在我想计算ABA，而不将Vector_A、Vector_B 转换回矩阵A, B。由于ABA 也是对称的，我想以与数组相同的方式保存ABA。我如何在 Julia 中做到这一点？

我还想计算X'AX，而不将Vector_A 转换回矩阵A，其中X' 表示transpose(X)。我如何在 Julia 中做到这一点？

【问题讨论】：

您确定需要高效存储矩阵吗？完整阵列的常规内存布局可能会超过性能优势。但是如果你真的需要节省内存，你最好写一个新的EfficientSymmetric{T} <: AbstractArray{T,2} 结构，这样你就可以“像数组一样使用它”。也许这个讨论是相关的？ discourse.julialang.org/t/symmetric-matrices/4086
AFAIK，Julia 具有稀疏矩阵类型。因此，与其构建自己的优化布局，不如使用稀疏矩阵。
我最后一次想到这一点是在我的ForecastEval.jl 包中处理算法模型置信度集时。我以两种方式实现该算法：1) 仅存储矩阵的上三角部分以节省内存，2) 存储整个矩阵，它使用两倍的内存，但允许我使用 BLAS 例程。第二种算法运行速度明显更快。第一个选项可能确实更适合您的问题，但至少也值得研究第二个。

标签： julia

【解决方案1】：

您需要实现自己的继承自AbstractMatrix 类型的数据结构。

例如，可以这样做：

struct SymmetricM{T} <: AbstractMatrix{T}
    data::Vector{T}
end

所以我们有一个对称矩阵，它只使用一个向量来存储它的数据。现在您需要实现函数，使其实际上表现得像一个矩阵，这样您就可以让 Julia 的魔法发挥作用。

我们首先提供新矩阵数据类型的大小。

function Base.size(m::SymmetricM) 
   n = ((8*length(m.data)+1)^0.5-1)/2
   nr = round(Int, n)
   @assert n ≈ nr "The vector length must match the number of triang matrix elements"
   (nr,nr)
end

在此代码中，nr 将每次在矩阵上完成 checkbounds 时计算。也许在您的生产实现中，您可能希望将其移动为SymmetricM 的字段。您会失去一些弹性并多存储 8 个字节，但会提高速度。

现在我们需要的下一个函数是根据矩阵索引计算向量的位置。这是一种可能的实现方式。

function getix(idx)::Int
    n = size(m)[1]
    row, col = idx
    #assume left/lower triangular
    if col > row
        row = col
        col = idx[1]
    end
    (row-1)*row/2 + col
end

现在我们可以实现getindex 和setindex 函数：

@inline function Base.getindex(m::SymmetricM, idx::Vararg{Int,2})
    @boundscheck checkbounds(m, idx...)
    m.data[getix(idx)]
end

@inline function Base.getindex(m::SymmetricM{T}, v::T, idx::Vararg{Int,2}) where T
    @boundscheck checkbounds(m, idx...)
    m.data[getix(idx)] = v
end

现在让我们测试一下：

julia> m = SymmetricM(collect(1:10))
4×4 SymmetricM{Int64}:
 1  2  4   7
 2  3  5   8
 4  5  6   9
 7  8  9  10

您可以看到我们只提供了一个三角形的元素（无论是下部还是上部 - 它们都是相同的） - 我们得到了完整的矩阵！

这确实是一个完全有效的 Julia 矩阵，所以所有的矩阵代数都应该适用于它：

julia> m * SymmetricM(collect(10:10:100))
4×4 Array{Int64,2}:
  700   840  1010  1290
  840  1020  1250  1630
 1010  1250  1580  2120
 1290  1630  2120  2940

请注意，乘法的结果是矩阵而不是 SymmetricM - 要获得 SymmetricM，您需要重载 * 运算符以接受 2 个 SymmetricM 参数。为了便于说明，让我们展示一个带有减号 - 的自定义运算符重载：

import Base.-
-(m1::SymmetricM, m2::SymmetricM) = SymmetricM(m1.data .- m2.data)

现在你会看到SymmetricM 的减法将返回另一个SymmetricM：

julia> m-m
4×4 SymmetricM{Int64}:
 0  0  0  0
 0  0  0  0
 0  0  0  0
 0  0  0  0

通过这种方式，您可以在 Julia 中构建完整的三角矩阵代数系统。

请注意，但是getix 函数有一个if 语句，因此在不使用data 字段的情况下访问SymmetricM 元素将比常规矩阵慢得多，因此也许您应该尝试重载尽可能多的矩阵您的项目所需的运算符。

【讨论】：