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
454 stars 122 forks source link

QR Decomposition works not fine #147

Closed KimKitsurag1 closed 5 months ago

KimKitsurag1 commented 5 months ago

Tried to use "arm_mat_qr_f32" and "arm_mat_qr_f64" to find eigenvalues by using QR method. Found out that R-matrix isn't triangular, because there are values that significantly differ from 0(~1e-1)

The code is following:

include "utils.h"

include

include "arm_math.h"

void printMatrix(float64_t matrix, uint32_t numRows, uint32_t numCols) { for (uint32_t i = 0; i < numRows; i++) { for (uint32_t j = 0; j < numCols; j++) { printf("%f\t", matrix[i numCols + j]); } printf("\n"); } } void eigen_values(arm_matrix_instance_f64 *A){ float64_t Q[9]={0}; float64_t R[9]={0}; arm_matrix_instance_f64 Q_mat = {.numCols = 3, .numRows = 3, .pData = Q}; arm_matrix_instance_f64 R_mat = {.numCols = 3, .numRows = 3, .pData = R};

 float64_t thresh = 10;
 float64_t tau[3]={0};
 float64_t A_[3]={0};
 float64_t B[3]={0};
 for (uint32_t k = 0; k < 100; k++){
     arm_mat_qr_f64(A, thresh, &R_mat, &Q_mat, tau, A_, B);
     printf("_______________________________");
     printf("\n");
     printMatrix(Q_mat.pData,3,3);
     printf("_______________________________");
     printf("\n");
     printMatrix(R_mat.pData,3,3);
     printf("_______________________________");
     printf("\n");
     arm_mat_mult_f64(&R_mat, &Q_mat, A);

     printMatrix(A->pData,3,3);
     printf("_______________________________");
     printf("\n");
 }

}

int main(){

 float64_t A_D[9] = {36.0f,   9.0f, -26.0f,   9.0f,  49.0f,  27.0f, -26.0f,  27.0f,  44.0f};
 arm_matrix_instance_f64 A = {.numCols = 3, .numRows =3, .pData = A_D};
 printMatrix(A_D, 3, 3);
 eigen_values(&A);
 printMatrix(A_D, 3, 3);

return 0;

} Снимок

christophe0606 commented 5 months ago

@KimKitsurag1 The matrix R is also containing the Householder reflectors below the diagonal like it is done in some matrix libraries. So indeed, there are some values below the diagonal and if you only want the R matrix part, you need to zero those values. (It is documented in description of the function).

For instance, first column below the diagonal in your screenshot is [0.110687, -0.319764] and it corresponds to Householder reflector [1,0.110687,-0.319764]

KimKitsurag1 commented 5 months ago

Thank you!