armadillo: matrix multiplication loss of precision

Question

I seem to experience a loss of precision when performing a matrix multiplication and would like to know how to prevent that. For example, assuming feat and beta are of the appropriate sizes,

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

The numbers I am working with are of rather smaller values, majority of numbers are less than 1e-3 so it is possible that what I am trying to achieve is impossible. I should also note that this is a MATLAB function calling a C++ function so the MEX compiler is involved. I do check the numbers in the mex function when they arrive and they are correct, it is only after the above line that I get incredibly wrong answers.

EDIT: I decide that it would not hurt to give full context of the program. Here is what I have so far. I have the line where the loss of precision comes up tagged with a comment.

EDIT2: Here are some examples of the matrices in question.
Feat_2d is 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 is 4609x4. I drop the last row for the multiplication of feat*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

Here is the code from 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
}

We need more info. Can you show some minimal examples of your matrices? also, how do you define "loss of precision"? How do you know you are actually loosing precision? If you are doing it properly, you can not compare it to what a computer does.... As maybe your comparison matrix is bad, not your own. — Ander Biguri
@AnderBiguri I made some edits showing some of the matrices, there seems to be a problem, in which Y results in basically zero operations. When I do the operation in matlab, the results are not the same. — John Halloran

Claes Rolen Claes Rolen · Accepted Answer · 2017-05-10T22:20:39

I can't see any direct errors in your code but I doubt that Armadillo has precision problems. I suspect that it might be something in your pointer to mat conversions. I would use the functions Armadillo provides in the .../mex_interface/armaMex.hpp file.

A simple example of a matrix multiplication (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;
}

and the .m file is

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)

giving the output

Matlab:
1.0e-06 *

0.240798243020273   0.410716970485213
0.559953800808707   0.974915983937571

Armadillo:
1.0e-06 *

0.240798243020273   0.410716970485213
0.559953800808707   0.974915983937571

armadillo: matrix multiplication loss of precision

1 Answers