Hello,
I am working on some D-optimal design widget that takes an input matrix, generates the square matrix, and calculates the determinant. This is what the determinant code looks like (omitting code that loads the matrices),

Code: Select all

#include <Eigen/Dense>
#include <Eigen/LU>

using namespace std;
using namespace Eigen;

// matrix for input PCA data
MatrixXd PCA_values, product_matrix;

// set row dimension to number of rows passed in call
current_number_of_rows = row_numbers.size();

// size matrices based on amount of data passed
PCA_values.resize(current_number_of_rows, num_data_columns);
product_matrix.resize(num_data_columns, num_data_columns);

// create the square matrix product
product_matrix = PCA_values.transpose() * PCA_values;

// calculate determinant
double current_determinant = 0.0;
current_determinant = product_matrix.determinant();


The answers I am getting don't make sense, so I am investigating if I need to do a transformation like FullPivLU. After reading some of the documentation, I modified the code to what is below in an attempt to perform a FullPivLU decomposition.

Code: Select all

#include <Eigen/Dense>
#include <Eigen/LU>

using namespace std;
using namespace Eigen;

// matrix for input PCA data
MatrixXd PCA_values, product_matrix;

// set row dimension to number of rows passed in call
current_number_of_rows = row_numbers.size();

// size matrices based on amount of data passed
PCA_values.resize(current_number_of_rows, num_data_columns);
product_matrix.resize(num_data_columns, num_data_columns);
FullPivLU<MatrixXd> lu(product_matrix);

// create the square matrix product
product_matrix = PCA_values.transpose() * PCA_values;

// calculate determinant
double current_determinant = 0.0;
current_determinant = lu.determinant();


The code compiles, but just runs for a while and segfaults. Can someone let me know if I have the syntax correct here? I can post the input file of the matrix for which the square is being calculated or the entire src, let me know.

LMHmedchem

Your first program looks correct to me. It would help if you could write some code which initializes PCA_values, ideally to a value for it is easy to compute the determinant but for which you get the wrong answer.

In your second program, you compute the LU decomposition of product_matrix before you initialize product_matrix. That will not work (though I would have expected that you got a nonsense answer instead of a segfault).

jitseniesen wrote:Your first program looks correct to me. It would help if you could write some code which initializes PCA_values, ideally to a value for it is easy to compute the determinant but for which you get the wrong answer.

In your second program, you compute the LU decomposition of product_matrix before you initialize product_matrix. That will not work (though I would have expected that you got a nonsense answer instead of a segfault).

Sorry about that, I just pasted the LU decomposition line in the wrong place. It does actually occur after the product matrix is generated. After trying a bit more, this will run with my untransformed input, but when I use a version of the same data that has been transformed to principle components, I get a segfault. I normally use PC because some of the columns are intercorrelated and the matrix can run sparse. This matrix is ~54% 0 cell values. That's not awfully sparse, but I tend to use PC just in case. The PCA version of the data has no intercorrelation and 0% 0 cell values. Is there any reason why I should get a segfault with the PCA input based on what I am doing below. This seems like pretty straightforward linear algebra.

LMHmedchem

Code: Select all

#include <Eigen/Dense>
#include <Eigen/LU>

using namespace std;
using namespace Eigen;

// matrix for input PCA data
MatrixXd PCA_values, product_matrix;

// set row dimension to number of rows passed in call
current_number_of_rows = row_numbers.size();

// size matrices based on amount of data passed
PCA_values.resize(current_number_of_rows, num_data_columns);
product_matrix.resize(num_data_columns, num_data_columns);

// create the square matrix product
product_matrix = PCA_values.transpose() * PCA_values;

//  take LU decomposition
FullPivLU<MatrixXd> lu(product_matrix);

// calculate determinant
double current_determinant = 0.0;
current_determinant = lu.determinant();


I don't see anything wrong but since this does not exactly correspond to your actual code, we cannot help. Please try to figure out the origin of the error by compiling in debug mode (-g, and do not define NDEBUG), and if needed, run your program within valgrind to spot memory errors.