如何在Julia中直接基于向量存储的对称矩阵计算ABA与X'AX

Let's break down the problem into two clear parts: calculating the scalar X'AX and computing the symmetric matrix ABA stored in the same upper-triangular column-major packed format as A and B.

Key Background: Index Mapping #

First, we need a helper function to map matrix indices (i,j) to their position in the packed vector. For an n×n symmetric matrix stored in upper-triangular column-major order:

• For any i ≤ j, the index is j*(j-1) ÷ 2 + i
• For i > j, we swap i and j (since the matrix is symmetric, A[i,j] = A[j,i])

function sym_index(i::Int, j::Int, n::Int)
    i > j && ((i, j) = (j, i))  # Swap to upper triangle if needed
    return j*(j-1) ÷ 2 + i
end

1. Calculating X'AX (Scalar) #

Since A is symmetric, we can compute the quadratic form directly using the packed vector without unpacking the full matrix. The formula leverages symmetry to avoid redundant calculations:

• Diagonal elements A[i,i] contribute A[i,i] * X[i]^2
• Off-diagonal elements A[i,j] (where i < j) contribute 2 * A[i,j] * X[i] * X[j] (since A[i,j] = A[j,i])

We iterate through the packed vector in its native column-major upper-triangular order for maximum efficiency:

function compute_XtAX(Vector_A::Vector{T}, X::Vector{T}, n::Int) where T<:Real
    result = zero(T)
    idx = 0
    for j in 1:n
        for i in 1:j
            idx += 1
            a_ij = Vector_A[idx]
            if i == j
                result += a_ij * X[i]^2
            else
                result += 2 * a_ij * X[i] * X[j]
            end
        end
    end
    return result
end

Example Verification #

For n=2, Vector_A = [1,2,3] (corresponding to A = [[1,2],[2,3]]), X = [7,8]:

compute_XtAX([1,2,3], [7,8], 2)  # Returns 465, matching manual calculation

2. Computing ABA (Packed Format) #

ABA is also a symmetric matrix, so we only need to compute its upper-triangular elements and store them in column-major order. We expand the matrix multiplication directly using the packed vectors, reducing time complexity to O(n³) (same as standard matrix multiplication):

function compute_ABA(Vector_A::Vector{T}, Vector_B::Vector{T}, n::Int) where T<:Real
    packed_len = n*(n+1) ÷ 2
    Vector_ABA = zeros(T, packed_len)
    idx_aba = 0
    
    for j in 1:n
        for i in 1:j
            idx_aba += 1
            val = zero(T)
            # Compute sum_{k} A[i,k] * (BA)[k,j]
            for k in 1:n
                ba_kj = zero(T)
                # Compute (BA)[k,j] = sum_{l} B[k,l] * A[l,j]
                for l in 1:n
                    b_kl = Vector_B[sym_index(k, l, n)]
                    a_lj = Vector_A[sym_index(l, j, n)]
                    ba_kj += b_kl * a_lj
                end
                a_ik = Vector_A[sym_index(i, k, n)]
                val += a_ik * ba_kj
            end
            Vector_ABA[idx_aba] = val
        end
    end
    return Vector_ABA
end

Example Verification #

For n=2, Vector_A = [1,2,3], Vector_B = [4,5,6] (corresponding to B = [[4,5],[5,6]]):

compute_ABA([1,2,3], [4,5,6], 2)  # Returns [48,79,130], matching manual calculation of ABA = [[48,79],[79,130]]

Notes #

• Efficiency: For large n, you can optimize further by adding multithreading to the inner k loop (using Julia's Threads.@threads) or precomputing frequently accessed elements.
• Generics: The functions use parametric types to support any real number type (e.g., Float32, Float64, Int).
• Memory Savings: Neither function reconstructs the full n×n matrices, making them memory-efficient for large datasets.

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