时间戳

首页 > 解决方案 > 如何计算将矩阵保存为向量的矩阵乘法

julia - 如何计算将矩阵保存为向量的矩阵乘法

问题描述

我有两个对称矩阵A和B一个向量X。的维数A是 n×n，的维数B是 n×n，的维数是 n×1。设矩阵的第 1 行第 1 列X的元素记为。ijAA[i,j]

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

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以与 matrix 相同的格式保存A。现在我想计算ABA而不转换Vector_A，Vector_B回到矩阵A, B。由于ABA也是对称的，我想以ABA与数组相同的方式保存。我怎样才能在 Julia 中做到这一点？

我还想在X'AX不转换Vector_A回表示的矩阵的A情况下进行计算。我怎样才能在 Julia 中做到这一点？X'transpose(X)

标签： julia

解决方案

您需要实现从该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

请注意，乘法的结果是 Matrix 而不是SymmetricM- 要获得 aSymmetricM您需要重载*运算符以接受 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语句，因此在不使用该字段的情况下访问SymmetricM元素data将比常规矩阵慢得多，因此也许您应该尝试重载项目所需的尽可能多的运算符。

推荐阅读