[eigen] [RFC] Test helper to make a matrix with exact rank.

[Keir Mierle <mierle@xxxxxxxxx>]

Here's a test routine to make a matrix with exact size and rank. It uses the SVD though, so may be too slow, and doesn't work on complex matrices (because SVD doesn't). Nevertheless, it is probably better than the current method of generating random matrices when it applies because the matrix will have excellent conditioning (i.e. the singular values will be 1, 1, 1, 1 ... (up to rank), 0, 0, 0, 0, 0).

It works by making a random (tall) matrix of the desired size, then setting the first 'rank' singular values to 1.0, and the remaining singular values to zero. Then the matrix is reconstituted by multiplying it back together, creating a matrix with the desired rows, columns, and rank.

Comments?

Keir
#include <Eigen/Core>
#include <Eigen/Cholesky>
#include <Eigen/LU>
#include <Eigen/QR>
#include <Eigen/SVD>
#include <Eigen/Array>

using namespace Eigen;
using namespace std;

// Make a well-conditioned matrix with precise rank. Must have rank <= cols <=
// rows. Takes a rows * cols SVD; this routine is not fast.
template<typename TMatrix>
void makeTallMatrix(int rows, int cols, int rank, TMatrix *matrix) {
  typedef typename TMatrix::Scalar Scalar;
  matrix->resize(rows, cols);
  matrix->setRandom();
  SVD<TMatrix> svd(*matrix);
  Matrix<Scalar, Dynamic, 1> s = svd.singularValues();
  int maxrank = std::min(rows, cols);
  assert(rank <= maxrank);
  for (int i = 0; i < rank; ++i)       s(i) = 1.0;
  for (int i = rank; i < maxrank; ++i) s(i) = 0.0;
  Matrix<Scalar, Dynamic, Dynamic> S(maxrank, maxrank);
  S.fill(0);
  S.block(0, 0, maxrank, maxrank).diagonal() = s;
  *matrix = svd.matrixU() * S * svd.matrixV().transpose();
}

template<typename TMatrix>
void makeMatrix(int rows, int cols, int rank, TMatrix *matrix) {
  if (rows < cols) {
    // Unfortunately a copy is necessary.
    TMatrix transposed = matrix->transpose();
    makeTallMatrix(cols, rows, rank, &transposed);
    *matrix = transposed.transpose();
  } else {
    makeTallMatrix(rows, cols, rank, matrix);
  }
}

template<typename TMatrix>
void makeMatrixWithRandomRank(TMatrix *matrix) {
  int rows = ei_random<int>(10,200);
  int cols = ei_random<int>(10,200);
  int rank = ei_random<int>(1, std::min(rows, cols));
  makeMatrix(rows, cols, rank, matrix);
}

int main() {
  MatrixXd m;
  makeMatrix(6, 4,  3, &m); cout << "m: \n" << m << endl;
  makeMatrix(23, 20,  6, &m); cout << "rank: " << m.lu().rank() << endl;
  makeMatrix(23, 20,  2, &m); cout << "rank: " << m.lu().rank() << endl;
  makeMatrix(20, 20, 11, &m); cout << "rank: " << m.lu().rank() << endl;
  makeMatrix(20, 23,  5, &m); cout << "rank: " << m.lu().rank() << endl;

  MatrixXf M;
  makeMatrix(23, 20,  6, &M); cout << "rank: " << M.lu().rank() << endl;
  makeMatrix(23, 20,  2, &M); cout << "rank: " << M.lu().rank() << endl;
  makeMatrix(20, 20, 11, &M); cout << "rank: " << M.lu().rank() << endl;
  makeMatrix(20, 23,  5, &M); cout << "rank: " << M.lu().rank() << endl;
  makeMatrixWithRandomRank(&M); cout << "rank: " << M.lu().rank() << " rows: " << M.rows() << " cols: " << M.cols() << endl;

  /*
   * Fails; SVD doesn't work for complex.
  MatrixXcf C;
  makeMatrix(23, 20,  6, &C); cout << "rank: " << C.lu().rank() << endl;
  makeMatrix(23, 20,  2, &C); cout << "rank: " << C.lu().rank() << endl;
  makeMatrix(20, 20, 11, &C); cout << "rank: " << C.lu().rank() << endl;
  makeMatrix(20, 23,  5, &C); cout << "rank: " << C.lu().rank() << endl;
  */
}