如何计算矩阵乘法，将矩阵保存为向量

Question

我有两个对称矩阵A，B和一个向量X。A的维数为n×n，B的维数为n×n，X的维数为n×1。设矩阵A的第i行、第j列的元素记为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中做到呢？

我还想在不将Vector_A转换回矩阵A的情况下计算X'AX，其中X'表示transpose(X)。我怎么才能在Julia中做到呢？

Answer 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

在此代码中，每次对矩阵执行checkbounds时，都会计算nr。也许在您的生产实现中，您可能希望将其移动到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元素将比访问常规矩阵慢得多，因此您可能应该尝试重载项目所需的尽可能多的运算符。