-1

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
rohlemax
  • 31
  • 4
  • I think you need more intermediate debuging displays. Not everyone can run your code and look at those values themselves. – hpaulj May 18 '22 at 23:45
  • ok I added some comments to make it clearer. I guess its more of a matrix-multiplication question than referring to the code. How can it generally be done; I just showed my first attempt, which fails as C != c – rohlemax May 19 '22 at 01:08
  • I vaguely recall from continuum mechanics studies a couple of decades ago that the 9 terms of a stress or strain matrix can be represented the 6 unique values. So I'm sure textbooks describe the matrix product operations on such vectors in detail, But don't expect a run-of-the-mill `numpy` programmer to be familiar with what you are doing. – hpaulj May 19 '22 at 01:59
  • ok ill put the question on mathematics; probably a better place for it. If done it on paper; it just breaks somewhere. Its a shame we cant use latex-style formulas here.. – rohlemax May 19 '22 at 12:43
  • You could explore it with `sympy` – hpaulj May 19 '22 at 14:32
  • I found my error. No I dont need sympy. the new matrix B can be constructed from P only. its equivalent, as the vector contains all the unique information of the symmetric matrix. – rohlemax May 19 '22 at 17:44

1 Answers1

0
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
rohlemax
  • 31
  • 4