ARM-software / CMSIS-DSP

CMSIS-DSP embedded compute library for Cortex-M and Cortex-A
https://arm-software.github.io/CMSIS-DSP
Apache License 2.0
564 stars 147 forks source link

Complex Matrix Linear Algebra #206

Open AvihooI opened 3 months ago

AvihooI commented 3 months ago

I apologize in advance, because this is not so much a CMSIS-DSP software issue (though it potentially could be), but rather a question about how to achieve certain functions/operations from linear algebra on complex matrices via CMSIS-DSP. If this is absolutely not the forum for this, then again, I apologize and I'd love to find what forum is more appropriate (I'm on the ARM Discord, though it's more generalist and the odds of me finding the help I need there are minimal).

I am dealing with a common DSP problem where I have complex fixed size square matrices that I need to perform a few common operations on:

  1. Multiplication comes right out of the box for complex matrices, which is appreciated.
  2. For finding the conjugate transpose (the Hermitian), I could re-interpret an NxN complex matrix as a 2Nx2N real matrix, transpose it and then multiply it by an appropriate diagonal matrix that negates the imaginary parts (this is effectively a diagonal matrix with 1s in odd places and -1s in even places). Is there a better way to achieve this? I noticed arm_mat_cmplx_trans_f32, but it's not obvious to me that it gives the Hermitian, or just the naive transpose.
  3. For inverting a square matrix (or applying any algebraic operation on it, really), there's a function that deals with real matrices, but none for complex matrices. Assuming my matrices are diagonalizable (which in my particular use case, they are), I could eigendecompose the matrix, invert all the eigenvalues and recompose using the eigenvector matrix and its Hermitian.
  4. Now this is where I'm really stuck: how do I achieve eigendecomposition? I could achieve it using QR decomposition, but there only exists a functionality for this for real matrices, not complex matrices.

What should I do on top of CMSIS-DSP to achieve that?

christophe0606 commented 3 months ago

@AvihooI Unfortunately, the complex matrix operations are quite limited in CMSIS-DSP and a lot are missing.

If a function is not available, it is generally not possible to make it from existing functions. It requires a new development.

There is an experimental branch : https://github.com/ARM-software/CMSIS-DSP/tree/cpp_complex

In this branch, we are adding a C++ API to make it easier to write complex matrix algorithms.

It is the C++ extension that is documented here : https://arm-software.github.io/CMSIS-DSP/latest/dsppp_main.html

The data type std::complex must be used.

Now, in C++ world, you could use Eigen library. I have had mixed results with performance on Cortex-M (because it looks like Eigen has been optimized more for bigger matrixes) but you'll find all the matrix algorithms you need with complex numbers but without any specific Cortex-M acceleration.

AvihooI commented 3 months ago

Oh, if I had to put my two cents on this, please don't couple this library with C++ 😟

Having said that, complex matrix eigendecomposition as I've recently discovered can be achieved by treating complex numbers as 2 by 2 real matrices that look like [a, -b],[b, a] for c = a + bi. And then matrix multiplication in this matrix space can achieve the same functionality as complex number multiplication.

As a consequence of this, N x N complex matrices can be represented as 2N x 2N matrices and so real eigendecomposition on that matrix will extract the the complex eigenvalues and eigenvectors as conjugate pairs.

This technique is obviously not memory and performance efficient, but it's a last resort solution for complex eigendecomposition if all you have are real matrix operations.

I have however decided to opt out of CMSIS-DSP and make my own implementation using NEON. If you want we can discuss these things further at some time. I do think a full fledged linear algebra library tailored specifically for ARM NEON has its place and I might take it upon myself to try implementing it.

christophe0606 commented 3 months ago

@AvihooI We have some Neon in CMSIS-DSP but obviously not enough and at some point we will add more.

C++ is optional in CMSIS-DSP.