[MLIR][Presburger] Add the related methods to support computing Matrix inverse

[kaitingwang]

Summary

Matrix inverse is useful when generating code from an affine schedule which must be full-rank (i.e. invertible).

[Groverkss]

Thank you for the patch!

Some comments on coding style to start off. I will put more comments on the algorithms used in a further review.

mlir/include/mlir/Analysis/Presburger/Matrix.h

48

48

Please mark this constructor as explicit.

185–196

Please use full stops at the end of sentences in comments.

mlir/lib/Analysis/Presburger/Matrix.cpp

25

vect[0] here can cause an out-of-bound access if vect is empty.

26

We prefer loops to be written with a pre-evaluated end condition and preincrement.

for (unsigned i = 0, e = vect.size(); i < e; ++i)

29–31

Braces can be dropped for single statement loops and if conditions:

for (unsigned j = 0; j < vect[0].size(); j++)
  at(i, j) = vect[i][j];

285

I don't think "quadratic" makes much sense. "square" seems to be a better term.

288–296

braces can be dropped for single statement if conditions and for loops.

300–310

Please use pre-increment for loops.

Also, the loop induction variables should be unsigned.

321–325

Braces can be dropped; use pre-evaluated for loop condition and pre-increment, please. Same for other loops as well.

331

quadratic -> square

363

This is better understood as "Matrix should be unimodular for integer inverse".

[kaitingwang]

@Groverkss Thanks for the detailed review. gtest does show the algorithm is tending on the expensive side. There was a reason that we went the route of computing the determinant/cofactor as suppose to a method such as LU decomposition. I can find out the precise reason from our scheduler person and update the MR tomorrow. Hope the style better matches your expectation and thanks for catching the zero-row vector bug.

[ RUN      ] MatrixTest.determinant
[       OK ] MatrixTest.determinant (106 ms)
[ RUN      ] MatrixTest.transpose
[       OK ] MatrixTest.transpose (0 ms)
[ RUN      ] MatrixTest.cofactor
[       OK ] MatrixTest.cofactor (107 ms)
[ RUN      ] MatrixTest.inverse
[       OK ] MatrixTest.inverse (111 ms)
[ RUN      ] MatrixTest.toVector
[       OK ] MatrixTest.toVector (0 ms)
[ RUN      ] MatrixTest.Matrix
[       OK ] MatrixTest.Matrix (108 ms)
[----------] 13 tests from MatrixTest (432 ms total)

[Groverkss]

In D133233#3773223, @kaitingwang wrote:

@Groverkss Thanks for the detailed review. gtest does show the algorithm is tending on the expensive side. There was a reason that we went the route of computing the determinant/cofactor as suppose to a method such as LU decomposition. I can find out the precise reason from our scheduler person and update the MR tomorrow. Hope the style better matches your expectation and thanks for catching the zero-row vector bug.

[ RUN      ] MatrixTest.determinant
[       OK ] MatrixTest.determinant (106 ms)
[ RUN      ] MatrixTest.transpose
[       OK ] MatrixTest.transpose (0 ms)
[ RUN      ] MatrixTest.cofactor
[       OK ] MatrixTest.cofactor (107 ms)
[ RUN      ] MatrixTest.inverse
[       OK ] MatrixTest.inverse (111 ms)
[ RUN      ] MatrixTest.toVector
[       OK ] MatrixTest.toVector (0 ms)
[ RUN      ] MatrixTest.Matrix
[       OK ] MatrixTest.Matrix (108 ms)
[----------] 13 tests from MatrixTest (432 ms total)

A simple reason I can see why the computation is expensive is that the determinant algorithm written is exponential, and each call makes a copy instead of creating a view. Fixing this would improve the performance quite a bit, but I would wait for the next round of review. Also, it would be helpful to know why you don't use LU decomposition (or Hermite normal form) before I can review the algorithm.

[kaitingwang]

[ RUN      ] MatrixTest.determinant
[       OK ] MatrixTest.determinant (106 ms)
[ RUN      ] MatrixTest.transpose
[       OK ] MatrixTest.transpose (0 ms)
[ RUN      ] MatrixTest.cofactor
[       OK ] MatrixTest.cofactor (107 ms)
[ RUN      ] MatrixTest.inverse
[       OK ] MatrixTest.inverse (111 ms)
[ RUN      ] MatrixTest.toVector
[       OK ] MatrixTest.toVector (0 ms)
[ RUN      ] MatrixTest.Matrix
[       OK ] MatrixTest.Matrix (108 ms)
[----------] 13 tests from MatrixTest (432 ms total)

A simple reason I can see why the computation is expensive is that the determinant algorithm written is exponential, and each call makes a copy instead of creating a view. Fixing this would improve the performance quite a bit, but I would wait for the next round of review. Also, it would be helpful to know why you don't use LU decomposition (or Hermite normal form) before I can review the algorithm.

We need to ensure/check that the determinant of the schedule matrix is always 1, -1 so that the inverse of the schedule can be non-fractional. This mostly stands from that we haven't figured out how to code generate for fractional schedule (which may require introducing if condition inside the newly generated loop to test for non-fraction, or use mod operator) We punted thinking about this to later. This leads us to pick (somewhat carelessly?) a method that implements finding out determinant.

mlir/lib/Analysis/Presburger/Matrix.cpp

25

Added a testcase for this.

[Groverkss]

From what I understand, you only need the integer inverse of a unimodular matrix. This is better accomplished by Hermite normal form.

The (column) hermite normal form for a matrix M is H = MU, such that U is unimodular, H is lower triangular, and the pivots in H are the max element in that column. Now, det(H) = 1 iff det(M) = 1, since det(U) = 1. det(H) = 1, actually forces H to be the identity matrix and U comes out to be inv(A). So you can simply find the hermite normal form, and use U as inv(A) if det(H) = 1 (since H is lower triangular this is just product of diagonals). Otherwise, M is not unimodular, and we cannot find an integer inverse.

We already have an algorithm to find the hermite normal form of a matrix upstream. I will refactor that algorithm to Matrix, and you can use that directly to get the inverse.

[Groverkss]

I sent the refactor here: https://reviews.llvm.org/D133510

In pseudo-code, what you would need to do to get the inverse is:

h, u = m.hermiteNormalForm()
check if all diagonal elements of h are 1 --> unimodularity check for matrix
if yes, return u
if no, return nothing since integer inverse does not exist

[Groverkss]

The patch was landed. You can use Hermite decomp to do inverse here.

[kaitingwang]

In D133233#3789860, @Groverkss wrote:

The patch was landed. You can use Hermite decomp to do inverse here.

Thanks. I'll close this pull request.

mlir/include/mlir/Analysis/Presburger/Matrix.h

public:
   /// The number of reserved rows and columns will be at least the number
   /// specified, and will always be sufficient to accomodate the number of rows
   /// and columns specified.
   ///
   /// Initially, the entries are initialized to ero.
   Matrix(unsigned rows, unsigned columns, unsigned reservedRows = 0,
          unsigned reservedColumns = 0);
 
+  /// Construct a matrix from a given vector of vector of integers.
+  explicit Matrix(const std::vector<std::vector<int64_t>> &vect);
+
   /// Return the identity matrix of the specified dimension.
   static Matrix identity(unsigned dimension);
 
   /// Access the element at the specified row and column.
   int64_t &at(unsigned row, unsigned column) {
     assert(row < nRows && "Row outside of range");
     assert(column < nColumns && "Column outside of range");
     return data[row * nReservedColumns + column];
public:
   /// Print the matrix.
   void print(raw_ostream &os) const;
   void dump() const;
 
   /// Return whether the Matrix is in a consistent state with all its
   /// invariants satisfied.
   bool hasConsistentState() const;
 
+  /// Return the determinant of the matrix.
+  int64_t determinant() const;
+  /// Return the transpose of the matrix.
+  Matrix transpose() const;
+  /// Return the cofactor of the matrix.
+  Matrix cofactor() const;
+  /// Return the inverse of the matrix, which must be square.
+  /// This method will not work for matrices without integer inverses.
+  Matrix inverse() const;
+
+  /// Return the Matrix as a vector of vector of integers
+  std::vector<std::vector<int64_t>> toVector() const;
+
 private:
   /// The current number of rows, columns, and reserved columns. The underlying
   /// data vector is viewed as an nRows x nReservedColumns matrix, of which the
   /// first nColumns columns are currently in use, and the remaining are
   /// reserved columns filled with zeros.
   unsigned nRows, nColumns, nReservedColumns;
 
   /// Stores the data. data.size() is equal to nRows * nReservedColumns.
   /// data.capacity() / nReservedColumns is the number of reserved rows.
   SmallVector<int64_t, 64> data;
 };
 
 } // namespace presburger
 } // namespace mlir
 
 #endif // MLIR_ANALYSIS_PRESBURGER_MATRIX_H

mlir/lib/Analysis/Presburger/Matrix.cpp

 Matrix::Matrix(unsigned rows, unsigned columns, unsigned reservedRows,
                unsigned reservedColumns)
     : nRows(rows), nColumns(columns),
       nReservedColumns(std::max(nColumns, reservedColumns)),
       data(nRows * nReservedColumns) {
   data.reserve(std::max(nRows, reservedRows) * nReservedColumns);
 }
 
+Matrix::Matrix(const std::vector<std::vector<int64_t>> &vect)
+    : Matrix::Matrix(vect.size(), vect.size() != 0 ? vect[0].size() : 0) {
+  for (unsigned i = 0, e = vect.size(); i < e; ++i) {
+    assert(vect[i].size() == nColumns &&
+           "2d array must be rectangular for conversion to Matrix.");
+    for (unsigned j = 0, e = vect[0].size(); j < e; ++j) {
+      at(i, j) = vect[i][j];
+    }
+  }
+}
+
 Matrix Matrix::identity(unsigned dimension) {
   Matrix matrix(dimension, dimension);
   for (unsigned i = 0; i < dimension; ++i)
     matrix(i, i) = 1;
   return matrix;
 }
 
 unsigned Matrix::getNumReservedRows() const {
bool Matrix::hasConsistentState() const {
   if (nColumns > nReservedColumns)
     return false;
   for (unsigned r = 0; r < nRows; ++r)
     for (unsigned c = nColumns; c < nReservedColumns; ++c)
       if (data[r * nReservedColumns + c] != 0)
         return false;
   return true;
 }
+
+int64_t Matrix::determinant() const {
+  assert(getNumRows() == getNumColumns() &&
+         "Matrix must be quadratic for determinant.");
+
+  unsigned dimension = getNumRows();
+  if (dimension == 0)
+    return 1;
+  if (dimension == 1)
+    return at(0, 0);
+  if (dimension == 2)
+    return at(0, 0) * at(1, 1) - at(0, 1) * at(1, 0);
+
+  int64_t result = 0;
+  int sign = 1;
+  for (unsigned i = 0; i < dimension; ++i) {
+    Matrix subVect(dimension - 1, dimension - 1);
+    for (unsigned m = 1; m < dimension; ++m) {
+      int z = 0;
+      for (unsigned n = 0; n < dimension; ++n) {
+        if (n != i) {
+          subVect.at(m - 1, z) = at(m, n);
+          z++;
+        }
+      }
+    }
+    result = result + sign * at(0, i) * subVect.determinant();
+    sign = -sign;
+  }
+  return result;
+}
+
+Matrix Matrix::transpose() const {
+  Matrix solution(getNumColumns(), getNumRows());
+  for (unsigned i = 0; i < nRows; ++i)
+    for (unsigned j = 0; j < nColumns; ++j)
+      solution.at(j, i) = at(i, j);
+  return solution;
+}
+
+Matrix Matrix::cofactor() const {
+  assert(getNumRows() == getNumColumns() &&
+         "Matrix must be quadratic for cofactor.");
+
+  Matrix solution(getNumRows(), getNumColumns());
+  Matrix subVect(getNumRows() - 1, getNumColumns() - 1);
+  for (unsigned i = 0; i < nRows; ++i) {
+    for (unsigned j = 0; j < nRows; ++j) {
+      int p = 0;
+      for (unsigned x = 0; x < nRows; ++x) {
+        if (x == i)
+          continue;
+        int q = 0;
+        for (unsigned y = 0; y < nRows; ++y) {
+          if (y == j)
+            continue;
+          subVect.at(p, q) = at(x, y);
+          q++;
+        }
+        p++;
+      }
+      solution.at(i, j) = pow(-1, i + j) * subVect.determinant();
+    }
+  }
+  return solution;
+}
+
+Matrix Matrix::inverse() const {
+  int64_t det = determinant();
+  assert((det == -1 || det == 1) && "Matrix must have integer inverse.");
+
+  Matrix solution(getNumRows(), getNumColumns());
+  solution = cofactor().transpose();
+  for (unsigned i = 0; i < nRows; ++i)
+    for (unsigned j = 0; j < nRows; ++j)
+      solution.at(i, j) = solution.at(i, j) * det;
+  return solution;
+}
+
+std::vector<std::vector<int64_t>> Matrix::toVector() const {
+  std::vector<std::vector<int64_t>> result(nRows,
+                                           std::vector<int64_t>(nColumns));
+  for (unsigned r = 0; r < nRows; ++r)
+    for (unsigned c = 0; c < nColumns; ++c)
+      result[r][c] = at(r, c);
+  return result;
+}

mlir/unittests/Analysis/Presburger/MatrixTest.cpp

TEST(MatrixTest, resize) {
   mat.resize(7, 7);
   ASSERT_TRUE(mat.hasConsistentState());
   EXPECT_EQ(mat.getNumRows(), 7u);
   EXPECT_EQ(mat.getNumColumns(), 7u);
   for (unsigned row = 0; row < 7; ++row)
     for (unsigned col = 0; col < 7; ++col)
       EXPECT_EQ(mat(row, col), row >= 3 || col >= 3 ? 0 : int(10 * row + col));
 }
+
+TEST(MatrixTest, determinant) {
+  Matrix mat1(3, 3);
+  std::vector<std::vector<int64_t>> test_data1(3, std::vector<int64_t>(3));
+  test_data1 = {{1, 99, 3}, {0, 1, 4}, {5, 6, 0}};
+  for (unsigned i = 0; i < 3; ++i)
+    for (unsigned j = 0; j < 3; ++j)
+      mat1.at(i, j) = test_data1[i][j];
+  EXPECT_EQ(mat1.determinant(), 1941);
+
+  Matrix mat2(4, 4);
+  std::vector<std::vector<int64_t>> test_data2(4, std::vector<int64_t>(4));
+  test_data2 = {{7, 3, 4, 1}, {0, 1, 1, 0}, {20, 6, 9, 2}, {6, 0, 1, 1}};
+  for (unsigned i = 0; i < 4; ++i)
+    for (unsigned j = 0; j < 4; ++j)
+      mat2.at(i, j) = test_data2[i][j];
+  EXPECT_EQ(mat2.determinant(), 1);
+
+  Matrix mat3(4, 4);
+  std::vector<std::vector<int64_t>> test_data3(4, std::vector<int64_t>(4));
+  test_data3 = {
+      {-1, -1, 8, 18}, {1, 2, 4, 4}, {10, 20, 41, 41}, {2, 2, -16, -35}};
+  for (unsigned i = 0; i < 4; ++i)
+    for (unsigned j = 0; j < 4; ++j)
+      mat3.at(i, j) = test_data3[i][j];
+  EXPECT_EQ(mat3.determinant(), -1);
+
+  Matrix mat_nonsquare(2, 1);
+  std::vector<std::vector<int64_t>> test_data_nonsquare(
+      2, std::vector<int64_t>(1));
+  test_data_nonsquare = {{1}, {2}};
+  for (unsigned i = 0; i < 2; ++i)
+    for (unsigned j = 0; j < 1; ++j)
+      mat_nonsquare.at(i, j) = test_data_nonsquare[i][j];
+  EXPECT_DEATH({ mat_nonsquare.determinant(); },
+               "Matrix must be quadratic for determinant.");
+}
+
+TEST(MatrixTest, transpose) {
+  Matrix mat(4, 2);
+  std::vector<std::vector<int64_t>> test_data(4, std::vector<int64_t>(2));
+  test_data = {{1, 2}, {3, 4}, {5, 6}, {7, 8}};
+  for (unsigned i = 0; i < 4; ++i)
+    for (unsigned j = 0; j < 2; ++j)
+      mat.at(i, j) = test_data[i][j];
+  std::vector<std::vector<int64_t>> test_res(2, std::vector<int64_t>(4));
+  test_res = {{1, 3, 5, 7}, {2, 4, 6, 8}};
+  for (unsigned i = 0; i < 2; ++i)
+    for (unsigned j = 0; j < 4; ++j)
+      EXPECT_EQ(mat.transpose().at(i, j), test_res[i][j]);
+}
+
+TEST(MatrixTest, cofactor) {
+  Matrix mat(4, 4);
+  std::vector<std::vector<int64_t>> test_data(4, std::vector<int64_t>(4));
+  test_data = {{7, 3, 4, 1}, {0, 1, 1, 0}, {20, 6, 9, 2}, {6, 0, 1, 1}};
+  for (unsigned i = 0; i < 4; ++i)
+    for (unsigned j = 0; j < 4; ++j)
+      mat.at(i, j) = test_data[i][j];
+  std::vector<std::vector<int64_t>> test_res(4, std::vector<int64_t>(4));
+  test_res = {{1, 8, -8, 2}, {-3, -17, 18, 0}, {0, -1, 1, -1}, {-1, -6, 6, 1}};
+  for (unsigned i = 0; i < 4; ++i)
+    for (unsigned j = 0; j < 4; ++j)
+      EXPECT_EQ(mat.cofactor().at(i, j), test_res[i][j]);
+
+  Matrix mat_nonsquare(2, 1);
+  std::vector<std::vector<int64_t>> test_data_nonsquare(
+      2, std::vector<int64_t>(1));
+  test_data_nonsquare = {{1}, {2}};
+  for (unsigned i = 0; i < 2; ++i)
+    for (unsigned j = 0; j < 1; ++j)
+      mat_nonsquare.at(i, j) = test_data_nonsquare[i][j];
+  EXPECT_DEATH({ mat_nonsquare.cofactor(); },
+               "Matrix must be quadratic for cofactor.");
+}
+
+TEST(MatrixTest, inverse) {
+  Matrix mat1(3, 3);
+  std::vector<std::vector<int64_t>> test_data1(3, std::vector<int64_t>(3));
+  test_data1 = {{1, 99, 3}, {0, 1, 4}, {5, 6, 0}};
+  for (unsigned i = 0; i < 3; ++i)
+    for (unsigned j = 0; j < 3; ++j)
+      mat1.at(i, j) = test_data1[i][j];
+  EXPECT_DEATH({ mat1.inverse(); }, "Matrix must have integer inverse.");
+
+  Matrix mat2(4, 4);
+  std::vector<std::vector<int64_t>> test_data2(4, std::vector<int64_t>(4));
+  test_data2 = {{7, 3, 4, 1}, {0, 1, 1, 0}, {20, 6, 9, 2}, {6, 0, 1, 1}};
+  for (unsigned i = 0; i < 4; ++i)
+    for (unsigned j = 0; j < 4; ++j)
+      mat2.at(i, j) = test_data2[i][j];
+  std::vector<std::vector<int64_t>> test_res2(4, std::vector<int64_t>(4));
+  test_res2 = {{1, -3, 0, -1}, {8, -17, -1, -6}, {-8, 18, 1, 6}, {2, 0, -1, 1}};
+  for (unsigned i = 0; i < 4; ++i)
+    for (unsigned j = 0; j < 4; ++j)
+      EXPECT_EQ(mat2.inverse().at(i, j), test_res2[i][j]);
+
+  Matrix mat3(4, 4);
+  std::vector<std::vector<int64_t>> test_data3(4, std::vector<int64_t>(4));
+  test_data3 = {
+      {-1, -1, 8, 18}, {1, 2, 4, 4}, {10, 20, 41, 41}, {2, 2, -16, -35}};
+  for (unsigned i = 0; i < 4; ++i)
+    for (unsigned j = 0; j < 4; ++j)
+      mat3.at(i, j) = test_data3[i][j];
+  std::vector<std::vector<int64_t>> test_res3(4, std::vector<int64_t>(4));
+  test_res3 = {
+      {38, -201, 20, 20}, {-19, 121, -12, -10}, {-2, -10, 1, -1}, {2, 0, 0, 1}};
+  for (unsigned i = 0; i < 4; ++i)
+    for (unsigned j = 0; j < 4; ++j)
+      EXPECT_EQ(mat3.inverse().at(i, j), test_res3[i][j]);
+}
+
+TEST(MatrixTest, toVector) {
+  Matrix mat(4, 4);
+  std::vector<std::vector<int64_t>> test_data(4, std::vector<int64_t>(4));
+  test_data = {{7, 3, 4, 1}, {0, 1, 1, 0}, {20, 6, 9, 2}, {6, 0, 1, 1}};
+  for (unsigned i = 0; i < 4; ++i)
+    for (unsigned j = 0; j < 4; ++j)
+      mat.at(i, j) = test_data[i][j];
+  ASSERT_TRUE(test_data == mat.toVector());
+}
+
+TEST(MatrixTest, Matrix) {
+  std::vector<std::vector<int64_t>> test_data1;
+  test_data1 = {{1, 2, 3, 4}, {1, 2, 3}, {1, 2, 3, 4}};
+  EXPECT_DEATH({ Matrix mat1(test_data1); },
+               "2d array must be rectangular for conversion to Matrix.");
+
+  std::vector<std::vector<int64_t>> test_data2(4, std::vector<int64_t>(4));
+  test_data2 = {{7, 3, 4, 1}, {0, 1, 1, 0}, {20, 6, 9, 2}, {6, 0, 1, 1}};
+  Matrix mat2(test_data2);
+  for (unsigned i = 0; i < 4; ++i)
+    for (unsigned j = 0; j < 4; ++j)
+      EXPECT_EQ(mat2.at(i, j), test_data2[i][j]);
+  std::vector<std::vector<int64_t>> test_data3; // empty vector.
+  Matrix mat3(test_data3);
+  EXPECT_EQ(mat3.getNumRows(), (unsigned)0);
+  EXPECT_EQ(mat3.getNumColumns(), (unsigned)0);
+}