[GiNaC-list] [patch] Improve automatic algorithm selection in matrix::solve()

[Vitaly Magerya]

Hi, folks. Since we've recently added solve_algo::markowitz, I
thought it would be a good idea to update the automatic elimination
algorithm selection logic we have for matrix::solve() (currently
living in matrix::echelon_form()).

To this end I've constructed a benchmark that generates random-ish
matrices and measures how fast matrix::solve() can solve them.
The matrix generation works by starting with a known simple
solution and then permuting and mixing rows and columns randomly.
You can find the code for this procedure at [1].

With that code I've generated a bunch of data points and run a
quick analysis on it. Overall, the fastest (on average) algorithm
seems to be:

    For numeric matrices:
        If ncells > 200 && density < 0.5: Markowitz
        Otherwise: Gauss

    For symbolic matrices:
        If density > 0.6 && ncells <= 12: Divfree
        If density > 0.6 && ncells < 120: Bareiss
        Otherwise: Markowitz

You can find the details (and some pretty charts) at [2].

The attached patch implements this recommendation.

Note that matrix::determinant() also has a similar automatic
algorithm selection logic, but I didn't touch it for now.

[1] https://github.com/magv/ginac-matrix-solve-bench/blob/master/mx-solve-bench.cpp
[2] https://github.com/magv/ginac-matrix-solve-bench/blob/master/analysis.ipynb
-------------- next part --------------
From 7c3c5284314b4c0edc4d23ca3ee4db11aa0839c9 Mon Sep 17 00:00:00 2001
From: Vitaly Magerya <magv at tx97.net>
Date: Thu, 14 Jun 2018 14:31:28 +0200
Subject: [PATCH] Improve automatic algorithm selection in matrix::solve()

---
 ginac/matrix.cpp | 35 +++++++++++++++++++++++++----------
 1 file changed, 25 insertions(+), 10 deletions(-)

diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp
index 3a01da2d..1bb33e0c 100644
--- a/ginac/matrix.cpp
+++ b/ginac/matrix.cpp
@@ -1211,21 +1211,36 @@ matrix::echelon_form(unsigned algo, int n)
 	if (algo == solve_algo::automatic) {
         // Gather some statistical information about the augmented matrix:
 		bool numeric_flag = true;
-		for (auto & r : m) {
+		unsigned ncells = col*row;
+		unsigned density = 0;
+		for (auto && r : m) {
 			if (!r.info(info_flags::numeric)) {
 				numeric_flag = false;
 				break;
 			}
+			density += !r.is_zero();
+		}
+		if (numeric_flag) {
+			// For numerical matrices Gauss is good, but Markowitz becomes
+			// better for large sparse matrices.
+			if ((ncells > 200) && (density < ncells/2)) {
+				algo = solve_algo::markowitz;
+			} else {
+				algo = solve_algo::gauss;
+			}
+		} else {
+			// For symbolic matrices Markowitz is good, but Bareiss/Divfree
+			// is better for small and dense matrices.
+			if ((ncells < 120) && (density*5 > ncells*3)) {
+				if (ncells <= 12) {
+					algo = solve_algo::divfree;
+				} else {
+					algo = solve_algo::bareiss;
+				}
+			} else {
+				algo = solve_algo::markowitz;
+			}
 		}
-		// Bareiss (fraction-free) elimination is generally a good guess:
-		algo = solve_algo::bareiss;
-		// For row<3, Bareiss elimination is equivalent to division free
-		// elimination but has more logistic overhead
-		if (row<3)
-			algo = solve_algo::divfree;
-		// This overrides any prior decisions.
-		if (numeric_flag)
-			algo = solve_algo::gauss;
 	}
 	// Eliminate the augmented matrix:
 	std::vector<unsigned> colid(col);
-- 
2.13.6