tatou: multiplication matricielle perte de précision

Question

Je semble éprouver une perte de précision lors de l'exécution d'une multiplication matricielle et aimerais savoir comment l'éviter. Par exemple, en supposant exploit et bêta sont des tailles appropriées,

Y = feat*beta.rows(0,N); 

Les chiffres avec lesquels je travaille sont des valeurs assez petites, la majorité des chiffres sont inférieurs à 1E-3 donc il est possible que ce que je suis en train accomplir est impossible. Je devrais également noter qu'il s'agit d'une fonction MATLAB appelant une fonction C++ de sorte que le compilateur MEX est impliqué. Je vérifie les chiffres dans la fonction mex quand ils arrivent et ils sont corrects, c'est seulement après la ligne ci-dessus que je reçois des réponses incroyablement fausses.

EDIT: Je décide que ça ne ferait pas de mal de donner le contexte complet du programme. Voici ce que j'ai jusqu'ici. J'ai la ligne où la perte de précision apparaît marquée d'un commentaire.

EDIT2: Voici quelques exemples des matrices en question.
Feat_2d est 5x4608

 0   0   0   0   0 
    0   0   0   0   0 
    0   0   0   0 4.8146 
    0   0 19.0266   0   0 
    0   0   0   0   0 

Beta_2d est 4609x4. Je laisse tomber la dernière ligne pour la multiplication de don * Beta_2d

-7.1486e-05 -1.6801e-04 1.0970e-05 3.7837e-04 
-8.7524e-05 1.8275e-04 -6.7857e-04 2.6267e-04 
-9.1812e-05 -6.5495e-05 -1.7687e-03 -3.2168e-04 
    0e+00  0e+00  0e+00  0e+00 
-4.5089e-04 -5.6013e-05 1.4841e-04 2.4912e-04 

Y =

6.8995e-310  0e+00 4.7430e-322 1.7802e-306 
6.8995e-310  0e+00 4.9407e-322 1.4463e-307 
0e+00  0e+00  0e+00 1.4463e-307 
0e+00  0e+00 1.2352e-322 1.2016e-306 
6.8996e-310 6.8996e-310 6.8995e-310 1.7802e-306 

Voici le code de EDIT1

#include <mex.h> 
#include <iostream> 
#include <armadillo> 

using namespace arma; 

void predict_bbox_reg(double *beta, int beta_dim[2], double *t_inv, int tinv_dim[2], double mu, double *feat, int feat_dim[2], double *boxes, int box_dim[2]){ 


    //convert pointers 
    //beta 
    arma::mat beta_2d = mat(beta_dim[0], beta_dim[1]); 
    for(int i = 0; i<beta_dim[0]; i++){ 
     for(int j = 0; j<beta_dim[1]; j++){ 
      beta_2d(i, j) = beta[i+beta_dim[0]*j]; 
     } 
    } 

    //t_inv 
    arma::mat tinv_2d = mat(tinv_dim[0], tinv_dim[1]); 
    for(int i = 0; i<tinv_dim[0]; i++){ 
     for(int j = 0; j<tinv_dim[1]; j++){ 
      tinv_2d(i, j) = t_inv[i+tinv_dim[0]*j]; 
     } 
    } 
    //feadoublet_2d 
    arma::mat feat_2d = mat(feat_dim[0], feat_dim[1]); 
    for(int i = 0; i<feat_dim[0]; i++){ 
     for(int j = 0; j<feat_dim[1]; j++){ 
      feat_2d(i, j) = feat[i+feat_dim[0]*j]; 
     } 
    } 

    //boxes 
    arma::mat box_2d = mat(box_dim[0], box_dim[1]); 
    for(int i = 0; i<box_dim[0]; i++){ 
     for(int j = 0; j<box_dim[1]; j++){ 
      box_2d(i, j) = boxes[i+box_dim[0]*j]; 
     } 
    } 


    arma::mat Y = mat(feat_dim[0], beta_dim[1]); 

    Y = feat_2d*beta_2d.rows(0,beta_dim[0]-2);// this is the precision loss 

    arma::mat y1 = beta_2d.row(beta_2d.n_rows-1); 
    Y.each_row() += y1.row(0); 

    //RETURNS SOMETHING 


}   

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { 

    int M = mxGetM(prhs[0]); 
    int N = mxGetN(prhs[0]); 

    int beta_dim[2] = {M,N}; 
    double *beta = mxGetPr(prhs[0]); 


    M = mxGetM(prhs[1]); 
    N = mxGetN(prhs[1]); 
    int tinv_dim[2] = {M,N}; 
    double *t_inv = mxGetPr(prhs[1]); 


    double mu = *mxGetPr(prhs[2]);  

    M = mxGetM(prhs[3]); 
    N = mxGetN(prhs[3]); 

    int feat_dim[2] = {M,N}; 
    double *feat = mxGetPr(prhs[3]); 

    M = mxGetM(prhs[4]); 
    N = mxGetN(prhs[4]); 

    int box_dim[2] = {M,N}; 
    double *ex_boxes = mxGetPr(prhs[4]); 

    predict_bbox_reg(beta, beta_dim, t_inv, tinv_dim, 
      mu, feat, feat_dim, ex_boxes, box_dim); 

    //RETURNS results to matlab 
} 

Answer 1

Je ne vois aucune erreur directe dans votre code mais je doute que Armadillo a des problèmes de précision. Je soupçonne que cela pourrait être quelque chose dans votre pointeur vers les conversions mat. J'utiliserais les fonctions fournies par Armadillo dans le fichier .../mex_interface/armaMex.hpp.

Un exemple simple d'une multiplication de matrice (mult_test.cpp):

#include <armadillo> 
#include "armaMex.hpp" 

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) 
{  
    // Convert to Armadillo 
    mat x1 = conv_to<mat>::from(armaGetPr(prhs[0],true)); 
    mat x2 = conv_to<mat>::from(armaGetPr(prhs[1],true)); 
    mat y(x1.n_rows,x2.n_cols); 

    // Do your stuff here: 
    y = x1*x2; 

    // Convert back to Matlab 
    plhs[0] = armaCreateMxMatrix(y.n_rows, y.n_cols, mxDOUBLE_CLASS, mxREAL); 
    armaSetPr(plhs[0], conv_to<mat>::from(y)); 
    return; 
} 

et le fichier .m est

X1 = [1 2 3 ; 4 5 6]; 
X2 = 1e-7*rand(3,2); 

Y = mult_test(X1,X2); 

disp('Matlab:') 
disp(X1*X2) 
disp('Armadillo:') 
disp(Y) 

donnant la sortie

Matlab: 
1.0e-06 * 

0.240798243020273 0.410716970485213 
0.559953800808707 0.974915983937571 

Armadillo: 
1.0e-06 * 

0.240798243020273 0.410716970485213 
0.559953800808707 0.974915983937571