构造向量化矩阵的等价变换

Question

矢量化解的等价变换

对于给定对称的4x4矩阵Q和3x4矩阵P，通过

得到3x3矩阵C。

C= p @ Q @ pt

可以证明，输出C将再次对称。同样的问题可以只用Q和C中唯一的元素利用它们的对称性来表述。为此，将矩阵矢量化，如下所示。

我想构造一个矩阵B，它将矢量化的矩阵彼此映射成这样:

c = B @ q

B必须是6x10，并且只能从P构造。如何从P得到B ?

我试过这个，但它似乎不工作。也许有人遇到过类似的问题?

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}")

输出:

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