《算法》BEYOND 部分程序 part 1

▶ 书中第六章部分程序，加上自己补充的代码，包括高斯消元法求解线性方程组，高斯 - 约旦消元法求解线性方程组

● 高斯消元法求解线性方程组，将原方程转化为上三角矩阵，然后从最后一个方程开始求解

 package package01;

 import edu.princeton.cs.algs4.StdOut;
 import edu.princeton.cs.algs4.StdRandom;

 public class class01
 {
     private static final double EPSILON = 1e-8;
     private final int m;                            // 方程个数
     private final int n;                            // 变量个数
     private double[][] a;                           // 增广矩阵
     private double[] x;                             // 方程的解
     private boolean haveSolution;                   // 方程是否有解

     public class01(double[][] A, double[] b)
     {
         m = A.length;
         n = A[0].length;
         if (b.length != m)
             throw new IllegalArgumentException("Dimensions disagree");
         a = new double[m][n + 1];
         x = new double[n];

         for (int i = 0; i < m; i++)                 // 搬运 A 和 b 到 a 中
         {
             for (int j = 0; j < n; j++)
                 a[i][j] = A[i][j];
         }
         for (int i = 0; i < m; i++)
             a[i][n] = b[i];

         for (int p = 0; p < Math.min(m, n); p++)    // 消元计算，a[p][p] 为当前主元
         {
             int max = p;                            // 从第 p 行往下寻找主元最大的行
             for (int i = p + 1; i < m; i++)
                 max = (Math.abs(a[i][p]) > Math.abs(a[max][p])) ? i : max;

             double[] temp = a[p];                   // 将最大主元行换到第 p 行
             a[p] = a[max];
             a[max] = temp;

             if (Math.abs(a[p][p]) <= EPSILON)       // 主元为 0，退化
                 continue;

             for (int i = p + 1; i < m; i++)         // 用主元消去后面的所有行
             {
                 double alpha = a[i][p] / a[p][p];
                 for (int j = p; j <= n; j++)
                     a[i][j] -= alpha * a[p][j];
             }
         }

         for (int i = Math.min(n - 1, m - 1); i >= 0; i--)   // 从最后一个方程开始求解 x[i]
         {
             double sum = 0.0;
             for (int j = i + 1; j < n; j++)
                 sum += a[i][j] * x[j];
             if (Math.abs(a[i][i]) > EPSILON)
                 x[i] = (a[i][n] - sum) / a[i][i];           // 解 x[i]
             else if (Math.abs(a[i][n] - sum) > EPSILON)     // x[i] 系数为 0，但是方程右边非 0，无解
             {
                 StdOut.println("No solution");
                 return;
             }
         }
         for (int i = n; i < m; i++)                         // m > n 的情形，检查剩余的方程是否有矛盾
         {
             double sum = 0.0;
             for (int j = 0; j < n; j++)
                 sum += a[i][j] * x[j];
             if (Math.abs(a[i][n] - sum) > EPSILON)
             {
                 StdOut.println("No solution");
                 return;
             }
         }

         for (int i = 0; i < m; i++)                         // 检验 A*x = b
         {
             double sum = 0.0;
             for (int j = 0; j < n; j++)
                 sum += A[i][j] * x[j];
             if (Math.abs(sum - b[i]) > EPSILON)
                 StdOut.println("not feasible\n b[" + i + "] = " + b[i] + ", sum = " + sum);
         }
         haveSolution = true;
     }

     public double[] solution()
     {
         return haveSolution ? x : null;
     }

     private static void test(String name, double[][] A, double[] b)
     {
         StdOut.println("----------------------------------------------------\n" + name + "\n----------------------------------------------------");
         class01 gaussian = new class01(A, b);
         double[] x = gaussian.solution();
         if (x != null)
         {
             for (int i = 0; i < x.length; i++)
                 StdOut.printf("%.6f\n", x[i]);
         }
         StdOut.println();
     }

     private static void test1()             // 3 × 3 非退化矩阵，x = [-1, 2, 2]
     {
         double[][] A = { { 0, 1,  1 },{ 2, 4, -2 },{ 0, 3, 15 } };
         double[] b = { 4, 2, 36 };
         test("test 1 (3-by-3 system, nonsingular)", A, b);
     }

     private static void test2()             // 3 × 3 无解
     {
         double[][] A = { { 1, -3,   1 },{ 2, -8,   8 },{ 3, -11, 9 } };
         double[] b = { 4, -2, 9 };
         test("test 2 (3-by-3 system, nonsingular)", A, b);
     }

     private static void test3()             // 5 × 5 无解
     {
         double[][] A = { { 2, -3, -1,  2,  3 },{ 4, -4, -1,  4, 11 },{ 2, -5, -2,  2, -1 },{ 0,  2,  1,  0,  4 },{ -4,  6,  0,  0,  7 } };
         double[] b = { 4, 4, 9, -6, 5 };
         test("test 3 (5-by-5 system, no solutions)", A, b);
     }

     private static void test4()             // 5 × 5 无穷多组解，x = [-6.25, -4.5, 0, 0, 1]
     {
         double[][] A = { { 2, -3, -1,  2,  3 },{ 4, -4, -1,  4, 11 },{ 2, -5, -2,  2, -1 },{ 0,  2,  1,  0,  4 },{ -4,  6,  0,  0,  7 } };
         double[] b = { 4, 4, 9, -5, 5 };
         test("test 4 (5-by-5 system, infinitely many solutions)", A, b);
     }

     private static void test5()             // 4 × 3 列满秩多解，x = [-1, 2, 2, ?]
     {
         double[][] A = { { 0, 1,  1 },{ 2, 4, -2 },{ 0, 3, 15 },{ 2, 8, 14 } };
         double[] b = { 4, 2, 36, 42 };
         test("test 7 (4-by-3 system, full rank)", A, b);
     }

     private static void test6()             // 4 × 3 列满秩无解
     {
         double[][] A = { { 0, 1,  1 },{ 2, 4, -2 },{ 0, 3, 15 },{ 2, 8, 14 } };
         double[] b = { 4, 2, 36, 40 };
         test("test 8 (4-by-3 system, no solution)", A, b);
     }

     private static void test7()             // 3×4 满秩，x = [3, -1, -2, 0]
     {
         double[][] A = { { 1, -3,   1,  1 },{ 2, -8,   8,  2 },{ -6,  3, -15,  3 } };
         double[] b = { 4, -2, 9 };
         test("test 9 (3-by-4 system, full rank)", A, b);
     }

     public static void main(String[] args)
     {
         test1();
         test2();
         test3();
         test4();
         test5();
         test6();
         test7();
         int n = Integer.parseInt(args[0]);  // 随机测试
         double[][] A = new double[n][n];
         for (int i = 0; i < n; i++)
         {
             for (int j = 0; j < n; j++)
                 A[i][j] = StdRandom.uniform(1000);
         }
         double[] b = new double[n];
         for (int i = 0; i < n; i++)
             b[i] = StdRandom.uniform(1000);
         test(n + "-by-" + n + " (probably nonsingular)", A, b);
     }
 }

● 高斯 - 约旦消元法求解线性方程组，将原方程转化为约旦标准型，无解时产生一个 yA == 0，yb != 0 的解

 package package01;

 import edu.princeton.cs.algs4.StdOut;
 import edu.princeton.cs.algs4.StdRandom;

 public class class01
 {
     private static final double EPSILON = 1e-8;
     private final int n;                        // 方程个数等于未知数个数
     private double[][] a;
     private double[] x;
     private boolean haveSolution = true;        // 方程是否有解

     public class01(double[][] A, double[] b)
     {
         n = b.length;
         a = new double[n][n + n + 1];
         x = new double[n];

         for (int i = 0; i < n; i++)             // 搬运
         {
             for (int j = 0; j < n; j++)
                 a[i][j] = A[i][j];
         }
         for (int i = 0; i < n; i++)             // 有解时 a[0:n-1][n:n*2] 最终为 A 的逆，否则某行为扩展解
             a[i][n + i] = 1.0;
         for (int i = 0; i < n; i++)
             a[i][n + n] = b[i];

         for (int p = 0; p < n; p++)                             // 消元计算
         {
             int max = p;
             for (int i = p + 1; i < n; i++)
                 max = (Math.abs(a[i][p]) > Math.abs(a[max][p])) ? i : max;

             double[] temp = a[p];
             a[p] = a[max];
             a[max] = temp;

             if (Math.abs(a[p][p]) <= EPSILON)
                 continue;

             for (int i = 0; i < n; i++)                         // 高斯-约旦消元
             {
                 double alpha = a[i][p] / a[p][p];
                 for (int j = 0; j <= n + n; j++)
                     a[i][j] -= (i != p) ? alpha * a[p][j] : 0;
             }
             for (int j = 0; j <= n + n; j++)                    // 主行归一化
                 a[p][j] /= (j != p) ? a[p][p] : 1;
             for (int i = 0; i < n; i++)                         // 被消去行的主元所在列元素强制归 0，主元归 1
                 a[i][p] = 0.0;
             a[p][p] = 1.0;
         }

         for (int i = 0; i < n; i++)                             // 求解 x
         {
             if (Math.abs(a[i][i]) > EPSILON)
                 x[i] = a[i][n + n] / a[i][i];
             else if (Math.abs(a[i][n + n]) > EPSILON)
             {
                 StdOut.println("No solution");
                 haveSolution = false;
                 break;
             }
         }
         if (!haveSolution)                              // 无解，输出约旦阵中第 i 行的结果
         {
             for (int i = 0; i < n; i++)
             {
                 if ((Math.abs(a[i][i]) <= EPSILON) && (Math.abs(a[i][n + n]) > EPSILON))
                 {
                     for (int j = 0; j < n; j++)
                         x[j] = a[i][n + j];
                 }
             }
         }

         if (haveSolution)                               // 检验结果
         {
             for (int i = 0; i < n; i++)                 // 有解则检验 A*x = b
             {
                 double sum = 0.0;
                 for (int j = 0; j < n; j++)
                     sum += A[i][j] * x[j];
                 if (Math.abs(sum - b[i]) > EPSILON)
                     StdOut.printf("not feasible\nb[%d] = %8.3f, sum = %8.3f\n", i, b[i], sum);
             }
         }
         else
         {
             for (int j = 0; j < n; j++)                 // 无解则检验 yA = 0, yb != 0
             {
                 double sum = 0.0;
                 for (int i = 0; i < n; i++)
                     sum += A[i][j] * x[i];
                 if (Math.abs(sum) > EPSILON)
                     StdOut.printf("invalid certificate of infeasibility\nsum = %8.3f\n", sum);
             }
             double sum = 0.0;
             for (int i = 0; i < n; i++)
                 sum += x[i] * b[i];
             if (Math.abs(sum) < EPSILON)
                 StdOut.printf("invalid certificate of infeasibility\nyb  = %8.3f\n", sum);
         }
     }

     public double[] solution()
     {
         return x;
     }

     public boolean haveActualSolution()
     {
         return haveSolution;
     }

     private static void test(String name, double[][] A, double[] b)
     {
         StdOut.println("----------------------------------------------------\n" + name + "\n----------------------------------------------------");
         class01 gaussian = new class01(A, b);
         double[] x = gaussian.solution();
         StdOut.println(gaussian.haveActualSolution() ? "Solution:" : "Certificate of infeasibility");
         for (int i = 0; i < x.length; i++)
             StdOut.printf("%10.6f\n", x[i]);
         StdOut.println();
     }

     private static void test1() // 3×3 非退化，x = [ -1, 2, 2 ]
     {
         double[][] A = { { 0, 1,  1 },{ 2, 4, -2 },{ 0, 3, 15 } };
         double[] b = { 4, 2, 36 };
         test("test 1", A, b);
     }

     private static void test2() // 3×3 无解，x = [ 1, 0, 1/3 ]
     {
         double[][] A = { { 2, -1,  1 },{ 3,  2, -4 },{ -6,  3, -3 } };
         double[] b = { 1, 4, 2 };
         test("test 5", A, b);
     }

     private static void test3() // 5×5 非退化，x = [ -6.25, -4.5, 0, 0, 1 ]
     {
         double[][] A = { { 2, -3, -1,  2,  3 },{ 4, -4, -1,  4, 11 },{ 2, -5, -2,  2, -1 },{ 0,  2,  1,  0,  4 },{ -4,  6,  0,  0,  7 } };
         double[] b = { 4, 4, 9, -5, 5 };
         test("test 4", A, b);
     }

     private static void test4() // 5×5 无解，x = [ -1, 0, 1, 1, 0 ]
     {
         double[][] A = { { 2, -3, -1,  2,  3 },{ 4, -4, -1,  4, 11 },{ 2, -5, -2,  2, -1 },{ 0,  2,  1,  0,  4 },{ -4,  6,  0,  0,  7 } };
         double[] b = { 4, 4, 9, -6, 5 };
         test("test 3", A, b);
     }

     private static void test5() // 3×3 无穷多组解，x = [ -1.375, 1.625, 0 ]
     {
         double[][] A = { { 1, -1,  2 },{ 4,  4, -2 },{ -2,  2, -4 } };
         double[] b = { -3, 1, 6 };
         test("test 6 (infinitely many solutions)", A, b);
     }

     public static void main(String[] args)
     {
         test1();
         test2();
         test3();
         test4();
         test5();
         int n = Integer.parseInt(args[0]);
         double[][] A = new double[n][n];
         double[] b = new double[n];
         for (int i = 0; i < n; i++)
         {
             for (int j = 0; j < n; j++)
                 A[i][j] = StdRandom.uniform(1000);
         }
         for (int i = 0; i < n; i++)
             b[i] = StdRandom.uniform(1000);
         test("random " + n + "-by-" + n + " (likely full rank)", A, b);

         for (int i = 0; i < n - 1; i++)
         {
             for (int j = 0; j < n; j++)
                 A[i][j] = StdRandom.uniform(1000);
         }
         for (int i = 0; i < n - 1; i++)
         {
             double alpha = StdRandom.uniform(11) - 5.0;
             for (int j = 0; j < n; j++)
                 A[n - 1][j] += alpha * A[i][j];
         }
         for (int i = 0; i < n; i++)
             b[i] = StdRandom.uniform(1000);
         test("random " + n + "-by-" + n + " (likely infeasible)", A, b);
     }
 }