Develop code and create and/or organize unit tests for matrix operations

[yungyuc]

This is to follow up issue #514. A set of matrix operations are listed in #502 (comment) . To facilitate the development, supply unit tests for the features that are already available, and add placeholders for those to be implemented.

Existing features or those may be easily extended from existing array operations:

• Matrix-matrix multiplication (naive) is implemented with PR Add GEMM implementation and tests for SimpleArray #587 Add a Python matmul/imatmul operator #596.
• Element-wise matrix multiplication.
  • Issue Implement element-wise matrix multiplication #595.
  • PR Implement element-wise addition and subtraction for SimpleArray #605.
• Matrix addition and subtraction (the same as element-wise addition).
  • Array-array PR Implemente Arithmetic operation for SimpleArray #533
  • Array-scalar PR Implement element-wise addition and subtraction for SimpleArray #605
• Element-wise matrix division.
  • Array-array PR Implemente Arithmetic operation for SimpleArray #533
  • Array-scalar PR Implement array to scalar division to SimpleArray #716
• Matrix transpose.
  • PRs Add the transpose function for SimpleArray #509 refactor function transpose and add the property .T #516.

The following items in the list will need to be developed:

• Generate identity matrix.
  • Issue Generate identify matrix #590.
  • PR Add identity matrix generation functionality to SimpleArray (#588) #591.
• Matrix power.
  • Issue Implement matrix power SimpleArray::pow() #821.
  • PR Add SimpleArray::pow() for matrix power #844.
• Matrix trace.
  • Issue Implement matrix trace SimpleArray::trace() #822.
  • PR Add member function SimpleArray::trace() #842.

Determinant is a complex task that should be tracked separately:

• Matrix determinant.
  • Issue Make the LU factorization helper class modmesh::detail::Lu stateful and outside detail #809 (LU decomposition) is a stepping stone.

The following items in the list will need underneath array implementation first:

• Element-wise equality test.
  • Issue Implement element-wise equality (eq) for SimpleArray #810.
  • PR Add element-wise eq() to SimpleArray #841.
• Scalar-matrix multiplication.

An item that is too big to be included in the scope of this issue:

• Matrix inversion: tracked in issue Implement a linear solver for SimpleArray #719.

[ThreeMonth03]

Hi, @5000user5000 , would you plan to implement inverse operation recently?
The algorithm of kalman filter need the inverse operation.

[yungyuc]

@ThreeMonth03 thank you for supplementing the dependency information. Could you please also supplement how you like the matrix inversion to be implemented? Additional information helps planning.

I am not certain about the plan of @5000user5000 . We do not need to assume any schedule for now.

[ThreeMonth03]

@ThreeMonth03 thank you for supplementing the dependency information. Could you please also supplement how you like the matrix inversion to be implemented? Additional information helps planning.

I am not certain about the plan of @5000user5000 . We do not need to assume any schedule for now.

The following alogrihtm is kalman filter:

Predict

$$x_{t|t-1} = F_t x_{t-1|t-1} + B_t u_t,$$ $$P_{t|t-1} = F_t P_{t-1|t-1} F_t^\top + Q_t.$$

Innovation

$$y_t = z_t - H_t x_{t|t-1},$$ $$S_t = H_t P_{t|t-1} H_t^\top + R_t.$$

Gain (via solve)

$$K_t = P_{t|t-1} H_t^\top S_t^{-1}.$$

Update

$$x_{t|t} = x_{t|t-1} + K_t y_t,$$ $$P_{t|t} = (I-K_t H_t) P_{t|t-1} (I-K_t H_t)^\top + K_t R_t K_t^\top.$$

To make the calculation precisely, the calculation of gain should be transformed like:

$$ S_t,K_t^\top ;=; \big(P_{t\mid t-1},H_t^\top\big)^\top. $$

Then we would solve the formula with LLT decomposition.

$$ K_t = \text{solve-with-LLT}\bigl(S_t,; (P_{t\mid t-1} H_t^\top)^\top\bigr)^\top. $$

[ThreeMonth03]

@ThreeMonth03 thank you for supplementing the dependency information. Could you please also supplement how you like the matrix inversion to be implemented? Additional information helps planning.
I am not certain about the plan of @5000user5000 . We do not need to assume any schedule for now.

The following alogrihtm is kalman filter:

Predict

x t | t − 1 = F t x t − 1 | t − 1 + B t u t , P t | t − 1 = F t P t − 1 | t − 1 F t ⊤ + Q t .

Innovation

y t = z t − H t x t | t − 1 , S t = H t P t | t − 1 H t ⊤ + R t .

Gain (via solve)

K t = P t | t − 1 H t ⊤ S t − 1 .

Update

x t | t = x t | t − 1 + K t y t , P t | t = ( I − K t H t ) P t | t − 1 ( I − K t H t ) ⊤ + K t R t K t ⊤ .

To make the calculation precisely, the calculation of gain should be transformed like:

S t , K t ⊤ ; = ; ( P t ∣ t − 1 , H t ⊤ ) ⊤ .

Then we would solve the formula with LLT decomposition.

K t = solve-with-LLT ( S t , ; ( P t ∣ t − 1 H t ⊤ ) ⊤ ) ⊤ .

I afterwards notice that the algorithm need the LLT solver instead of the inverse solver, it is possible for solving the aforementioned formula without solving inverse matrix.

[tkoyama010]

I will do Generate identity matrix.

[yungyuc]

I will do Generate identity matrix.

@tkoyama010 Issue #590 is created for you.