python - 构造向量化矩阵的等效变换

矢量化解的等价变换

对于给定的对称 4x4 矩阵 Q3x4 矩阵 P3x3 矩阵 C 通过

C=P@Q@P.T

可以证明输出 C 将再次对称。仅使用 QC 中的独特元素利用它们的对称性,可以制定相同的问题。为此,矩阵被向量化,如下所示。

我想构造一个矩阵 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

回答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

相似文章

r - 矩阵列组合算法

我尝试按R的矩阵列进行所有组合,如下所示首先,我有一个4*3矩阵,如下所示>a[,1][,2][,3][,4][1,]1011[2,]1000[3,]1111然后我们要列出总共312*2=12组合,例...

最新文章