Construct equivalent transform for vectorized Matrix

Question

Equivalent transform for vectorized solution

For a given symmetric 4x4 matrix Q and a 3x4 matrix P the 3x3 matrix C is obtained through

C=P @ Q @ P.T

It can be shown that the output C will be symmetric again. The same problem can be formulated using only the unique elements in Q and C exploiting their symmetry. To do so, the matrices are vectorized as seen below.

I want to construct a matrix B that maps the vectorized matrices onto each other like so:

c = B @ q

B must be a 6x10 and should be constructable from P only. How can I get B from P?

I tried this, but it doesnt seem to work. Maybe someone has experienced a similar problem?

import numpy as np


def vectorize(A, ord='c'):
    """
    Symmetric matrix to vector e.g:
    [[1, 2, 3],
     [2, 4, 5],
     [3, 5, 6]] -> [1, 2, 3, 4, 5, 6] (c-order, row-col)
                -> [1, 2, 4, 3, 5, 6] (f-order, col-row)
    """
    # upper triangle mask
    m = np.triu(np.ones_like(A, dtype=bool)).flatten(order=ord)
    return A.flatten(order=ord)[m]


def B(P):
    B = np.zeros((6, 10))
    counter = 0
    # the i,j entry in C depends on the i, j columns in P
    for i in range(3):
        for j in range(i, 3):
            coeffs = np.outer(P[i], P[j])
            B[counter] = vectorize(coeffs)
            counter += 1
    return B


if __name__ == '__main__':
    # original transform
    P = np.arange(12).reshape((3, 4))

    # calculated transform for vectorized matrix
    _B = B(P)

    # some random symmetric matrix
    Q = np.array([[1, 2, 3, 4],
                  [2, 5, 6, 7],
                  [3, 6, 8, 9],
                  [4, 7, 9, 10]])

    # if B is an equivilant transform to P, these should be similar
    C = P @ Q @ P.T
    c = _B @ vectorize(Q)
    print(f"q: {vectorize(Q)}\n"
          f"C: {vectorize(C)}\n"
          f"c: {c}")

Output:

q: [ 1  2  3  4  5  6  7  8  9 10]
C: [ 301  949 2973 1597 4997 8397]
c: [ 214  542  870 1946 3154 5438] <-- not the same

Answer 1

import numpy as np


def vec_from_mat(A, order='c'):
    """
    packs the unique elements of symmetric matrix A into a vector
    :param A: symmetric matrix
    :return:
    """
    return A[np.triu_indices(A.shape[0])].flatten(order=order)

def B_from_P(P):
    """
    returns a 6x10 matrix that maps the 10 unique elements of a symmetric 4x4 matrix Q on the 6 unique elements of a
    3x3 matrix C to linearize the equation C=PTQP to c=Bv
    :param P: 3x4 matrix
    :return: B with shape (6, 10)
    """
    n, m = P.shape
    b1, b2 = (n * (n + 1) // 2), (m * (m + 1) // 2)
    B = np.zeros((b1, b2))
    for a, (i, j) in enumerate(zip(*np.triu_indices(n))):
        coeffs = np.outer(P[i], P[j])
        # collect coefficients from lower and upper triangle of symmetric matrix
        B[a] = vec_from_mat(coeffs) + vec_from_mat(np.triu(coeffs.T, k=1))
    return B