关于使用克莱姆法则和矩阵消元法解线性方程组的效率的比较

《众所周知》，在线性代数这门课中，我们学习了两种解线性方程组的方法：克莱姆法则（Cramer's Rule）和矩阵消元法。但是在考研数学李永乐的课堂上了解到，由于计算量太大，克莱姆法则并不适合解线性方程组（亏我还在洪老师课堂上给学弟学妹们安利这种方法解方程组，我错了）。那么，使用克莱姆法则的效率到底能低到什么程度呢？本文将使用C#语言在Debug模式下分别使用两种方法解10000个含有30个未知数，30个方程的非齐次线性方程组，以比较它们的优劣。为了简化测评，本程序的测试数据均排除了方程组无解的情况，确保方程组有且仅有唯一解。

C#代码如下：

using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Text;

namespace ConsoleApp1
{
    class Program
    {
        class Dets
        {
            public double[,] xishu, zengguang;
            public double[][,] replace;
        }
        const int N = 30;
        static void Main(string[] args)
        {
            Random rand = new Random();
            List<Dets> dets = new List<Dets>();
            while (dets.Count < 10000)
            {
                Dets d = new Dets();
                d.xishu = new double[N, N];
                d.zengguang = new double[N, N + 1];
                for (int i = 0; i < N; i++)
                {
                    for (int j = 0; j < N; j++)
                    {
                        d.zengguang[i, j] = d.xishu[i, j] = rand.NextDouble();
                    }
                    d.zengguang[i, N] = rand.NextDouble();
                }
                if (getDet(d.xishu) == 0) continue;
                d.replace = new double[N][,];
                for (int k = 0; k < N; k++)
                {
                    d.replace[k] = new double[N, N];
                    for (int i = 0; i < N; i++)
                    {
                        for (int j = 0; j < N; j++)
                        {
                            d.replace[k][i, j] = d.xishu[i, j];
                        }
                    }
                    for (int i = 0; i < N; i++)
                    {
                        d.replace[k][i, k] = d.zengguang[i, N];
                    }
                }
                dets.Add(d);
            }
            StringBuilder sb = new StringBuilder();
            Stopwatch sw_klm = new Stopwatch(), sw_hzj = new Stopwatch();
            for (int k = 0; k < dets.Count; k++)
            {
                Console.WriteLine($"第{k + 1}组非齐次线性方程组=======================");
                sw_klm.Start();
                double d = getDet(dets[k].xishu, false);
                for (int i = 0; i < N; i++)
                {
                    double di = getDet(dets[k].replace[i], false);
                    sb.Append($"x{i + 1} = {di / d}\r\n");
                }
                sw_klm.Stop();
                Console.WriteLine("使用克莱姆法则答案：\r\n" + sb);
                sb.Clear();
                //=================================================================
                sw_hzj.Start();
                hangzuijian(dets[k].zengguang);
                for (int i = 0; i < N; i++)
                {
                    sb.Append($"x{i + 1} = {dets[k].zengguang[i, N]}\r\n");
                }
                sw_hzj.Stop();
                Console.WriteLine("使用初等行变换化行最简答案：\r\n" + sb);
                sb.Clear();
            }
            Console.WriteLine("使用克莱姆法则用时：" + sw_klm.ElapsedMilliseconds + "ms，使用初等行变换化行最简用时：" + sw_hzj.ElapsedMilliseconds + "ms");
        }
        static void printArray(double[,] a)
        {
            int row = a.GetLength(0), col = a.GetLength(1);
            for (int i = 0; i < row; i++)
            {
                for (int j = 0; j < col; j++)
                {
                    Console.Write(a[i, j] + " ");
                }
                Console.WriteLine();
            }
        }
        static void swap<T>(ref T x, ref T y)
        {
            T tmp = x;
            x = y;
            y = tmp;
        }
        static double getDet(double[,] mat, bool copy = true)
        {
            if (copy) mat = mat.Clone() as double[,];
            int n = mat.GetLength(0);
            double ans = 1;
            for (int i = 0; i < n; i++)
            {
                int j = i;
                while (j < n && mat[j, i] == 0) j++;
                if (j >= n) return 0;
                if (j > i)
                {
                    for (int k = i; k < n; k++)
                    {
                        swap(ref mat[i, k], ref mat[j, k]);
                    }
                }
                ans *= mat[i, i];
                //化为上三角（下三角归零）
                for (j = i + 1; j < n; j++)
                {
                    if (mat[j, i] == 0) continue;
                    for (int k = i + 1; k < n; k++)
                    {
                        mat[j, k] -= mat[j, i] / mat[i, i] * mat[i, k];
                    }
                    mat[j, i] = 0;
                }
            }
            return ans;
        }
        static void hangzuijian(double[,] mat)
        {
            int row = mat.GetLength(0), col = mat.GetLength(1);
            for (int i = 0; i < row; i++)
            {
                int j = i, startCol = i;
                while (startCol < col)
                {
                    while (j < row && mat[j, startCol] == 0) j++;
                    if (j >= row)
                    {
                        startCol++;
                        j = i;
                    }
                    else
                    {
                        break;
                    }
                }
                if (startCol >= col) break;
                if (j > i)
                {
                    for (int k = startCol; k < col; k++)
                    {
                        swap(ref mat[i, k], ref mat[j, k]);
                    }
                }
                //归一
                if (mat[i, startCol] != 1)
                {
                    for (j = startCol + 1; j < col; j++)
                    {
                        mat[i, j] /= mat[i, startCol];
                    }
                    mat[i, startCol] = 1;
                }
                //排他
                for (j = 0; j < row; j++)
                {
                    if (j == i || mat[j, startCol] == 0) continue;
                    for (int k = startCol + 1; k < col; k++)
                    {
                        mat[j, k] -= mat[j, startCol] * mat[i, k];
                    }
                    mat[j, startCol] = 0;
                }
            }
        }
    }
}

运行结果如图所示：

从结果图上可知，矩阵消元法的效率远高于克莱姆法则，但这仅限大量频繁解未知数较多的方程组。对于一些简单的方程组（例如二元二次方程组），前者则稍慢于后者。