CODEWITHSUNDEEP
julia
Mizzle
Mizzle

Reputation: 757

How to calculate matrix multiplication in which matrix is saved as vector

I have two symmetric matrices A, B and a vector X. The dimension of A is n by n, the dimension of B is n by n, the dimension of X is n by 1. Let the element at ith row and jth column of matrix A denoted by A[i,j].

Since A is symmetric, only each column of the upper triangular matrix of A is saved. The matrix A is saved as an array：

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]]

The matrix B is saved in the same format as matrix A. Now I would like to calculate ABA without transforming Vector_A, Vector_B back to matrix A, B. Since ABA is also symmetric, I would like to save the ABA in the same way as an array. How can I do it in Julia?

I would also like to calculate X'AX without transforming Vector_A back to matrix A where X' denotes transpose(X). How can I do it in Julia?

Upvotes: 3

Views: 444

Answers (1)

Przemyslaw Szufel
Przemyslaw Szufel

Reputation: 42214

You need to implement your own data structures that inherit from the the AbstractMatrix type.

For example this could be done as:

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

So we have a symmetric matrix that is using only a vector for its data storage. Now you need to implement functions so it actually behaves like a matrix so you can let the Julia magic work.

We start by providing the size of our new matrix datatype.

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

In this code nr will be calculate every time to checkbounds is done on matrix. Perhaps in your production implementation you might want to move it to be a field of SymmetricM. You would scarify some elasticity and store 8 bytes more but would gain on the speed.

Now the next function we need is to calculate position of the vector on the base of matrix indices. Here is one possible implementation.

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

Having that now we can implement getindex and setindex functions:

@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

Now let us test this thing:

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

You can see that we have provide elements of only one triangle (be it the lower or upper - they are the same) - and we got the full matrix!

This is indeed a fully valid Julia matrix so all matrix algebra should work on it:

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

Note that the result of multiplication is a Matrix rather than SymmetricM - to get a SymmetricM you need to overload the * operator to accept 2 SymmetricM arguments. For illustrative purposes let us show a custom operator overloading with the minus sign -:

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

And now you will see that substraction of SymmetricM is going to return another SymmetricM:

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

In this way you can build a full triangular matrix algebra system in Julia.

Note that however the getix function has an if statement so access to SymmetricM elements without using the data field will be much slower than those of a regular matrix so perhaps you should try to overload as many operators as is required for your project.

Upvotes: 1

Related Questions