C# 实现基本线性代数概念

引言

我在这里分享我的 C# 实现的基本线性代数概念。完整的代码已发布到 GitHub (https://github.com/elsheimy/Elsheimy.Components.Linears)，采用 MIT 许可证和 Nuget (https://nuget.net.cn/packages/Elsheimy.Components.Linears/)。您可以自由地修改和使用代码，没有任何限制或约束（同样，不提供任何形式的保证）。请让我知道您的反馈、意见、建议和更正。

代码涵盖以下内容

矩阵

• 矩阵加法/减法
• 矩阵连接/收缩/提取
• 矩阵化简和消元（高斯和高斯-若尔当）
• 矩阵秩、解和解的状态（唯一解、无穷解和无解）
• 使用拉普拉斯展开法计算行列式
• 可逆性和逆矩阵
• 矩阵镜像
• 矩阵乘法和幂
• 转置
• 单位矩阵
• 变换矩阵
• 投影矩阵
• 反射矩阵
• 2D 和 3D 旋转矩阵
• 2D 和 3D 错切矩阵
• 克隆和相等比较器

向量

• 向量加法/减法
• 向量之间的夹角
• 向量归一化和模长计算器
• 点积和叉积
• 向量缩放
• 投影到其他向量和子空间
• 克隆和相等比较器

该设计基于对 Array 对象进行操作的原始函数。为矩阵和向量（分别是 Matrix 和 Vector）创建了包装类，以隐藏这些核心函数。添加了静态初始化器和运算符，以简化包装类的使用。

引用

我不会讨论一个函数在数学上是如何工作的。尽管如此，欢迎您随时纠正我。

让我们简单地了解一下代码的一些亮点。

矩阵

矩阵消元

此函数演示了高斯消元法和高斯-若尔当消元法。它可以用于将矩阵简化为行阶梯形 (REF / Gaussian) 或简化行阶梯形 (RREF / Gauss-Jordan)。它可以用于求解具有多个增广列的矩阵，并计算矩阵的秩。

/// <summary>
/// Reduces matrix to row-echelon (REF/Gauss) or 
/// reduced row-echelon (RREF/Gauss-Jordan) form and solves for augmented columns (if any.)
/// </summary>
public static MatrixEliminationResult Eliminate(double[,] input, 
              MatrixReductionForm form, int augmentedCols = 0) {
  int totalRowCount = input.GetLength(0);
  int totalColCount = input.GetLength(1);

  if (augmentedCols >= totalColCount)
    throw new ArgumentException(nameof(augmentedCols),
      Properties.Resources.Exception_TooMuchAugmentedColumns);

  MatrixEliminationResult result = new MatrixEliminationResult();
  double[,] output = CreateCopy(input);
 
  // number of pivots found
  int numPivots = 0;
 
  // loop through columns, exclude augmented columns
  for (int col = 0; col < totalColCount - augmentedCols; col++) {
    int? pivotRow = FindPivot(output, numPivots, col, totalRowCount);
    if (pivotRow == null)
      continue; // no pivots, go to another column
    ReduceRow(output, pivotRow.Value, col, totalColCount);
    SwitchRows(output, pivotRow.Value, numPivots, totalColCount);
    pivotRow = numPivots;
    numPivots++;

    // Eliminate Previous Rows
    if (form == MatrixReductionForm.ReducedRowEchelonForm) {
      for (int tmpRow = 0; tmpRow < pivotRow; tmpRow++)
        EliminateRow(output, tmpRow, pivotRow.Value, col, totalColCount);
    }

    // Eliminate Next Rows
    for (int tmpRow = pivotRow.Value; tmpRow < totalRowCount; tmpRow++)
      EliminateRow(output, tmpRow, pivotRow.Value, col, totalColCount);
  }

  result.Rank = numPivots;
  result.FullMatrix = output;
  result.AugmentedColumnCount = augmentedCols;
  if (augmentedCols > 0 && 
      form == MatrixReductionForm.ReducedRowEchelonForm) { // matrix has solution
    result.Solution = ExtractColumns(output, totalColCount - augmentedCols, totalColCount - 1);
  }

  return result;
}

private static int? FindPivot(double[,] input, int startRow, int col, int rowCount) {
  for (int i = startRow; i < rowCount; i++) {
    if (input[i, col] != 0)
      return i;
  }

  return null;
}

private static void SwitchRows(double[,] input, int row1, int row2, int colCount) {
  if (row1 == row2)
    return;

  for (int col = 0; col < colCount; col++) {
    var tmpVal1 = input[row1, col];
    input[row1, col] = input[row2, col];
    input[row2, col] = tmpVal1;
  }
}

private static void ReduceRow(double[,] input, int row, int col, int colCount) {
  var coeffecient = 1.0 / input[row, col];

  if (coeffecient == 1)
    return;

  for (; col < colCount; col++)
    input[row, col] *= coeffecient;
}

/// <summary>
/// Eliminates row using another pivot row.
/// </summary>
private static void EliminateRow(double[,] input, int row, int pivotRow, 
                                 int pivotCol, int colCount) {
  if (pivotRow == row)
    return;

  if (input[row, pivotCol] == 0)
    return;

  double coeffecient = input[row, pivotCol];

  for (int col = pivotCol; col < colCount; col++) {
    input[row, col] -= input[pivotRow, col] * coeffecient;
  }
}

此函数已使用以下函数包装到我们的 Matrix 包装器中

public virtual Matrix Reduce(MatrixReductionForm form, int? augmentedColCount = null);
public virtual MatrixSolutionState GetSolutionState(int? augmentedColCount = null);
public virtual Matrix Solve(int? augmentedColCount = null);
public virtual int GetRank(int? augmentedColCount = null);

矩阵行列式

此函数行列式使用拉普拉斯展开法。它通过计算matrix的叉积来工作。这是实现

/// <summary>
/// Calculates determinant. Internally uses Laplace Expansion method.
/// </summary>
/// <remarks>
/// Returns 1 for an empty matrix. 
/// See https://math.stackexchange.com/questions/1762537/what-is-the-determinant-of/1762542
/// </remarks>
public static double Determinant(double[,] input) {
  var results = CrossProduct(input);

  return results.Sum();
}

public static double[] CrossProduct(double[,] input) {
  int rowCount = input.GetLength(0);
  int colCount = input.GetLength(1);

  if (rowCount != colCount)
    throw new InvalidOperationException(Properties.Resources.Exception_RequiredSquareMatrix);

  if (rowCount == 0)
    return new double[] { 1 };

  if (rowCount == 1)
    return new double[] { input[0, 0] };

  if (rowCount == 2)
    return new double[] { (input[0, 0] * input[1, 1]) - (input[0, 1] * input[1, 0]) };

  double[] results = new double[colCount];

  for (int col = 0; col < colCount; col++) {
    // checkerboard pattern, even col  = 1, odd col = -1
    var checkerboardFactor = col % 2.0 == 0 ? 1 : -1;
    var coeffecient = input[0, col];

    var crossMatrix = GetCrossMatrix(input, 0, col);
    results[col] = checkerboardFactor * (coeffecient * Determinant(crossMatrix));
  }

  return results;
}

/// <summary>
/// Retrieves all matrix entries except the specified row and col. 
/// Used in cross product and determinant functions.
/// </summary>
public static double[,] GetCrossMatrix(double[,] input, int skipRow, int skipCol) {
  int rowCount = input.GetLength(0);
  int colCount = input.GetLength(1);

  var output = new double[rowCount - 1, colCount - 1];
  int outputRow = 0;

  for (int row = 0; row < rowCount; row++) {
    if (row == skipRow)
      continue;

    int outputCol = 0;

    for (int col = 0; col < colCount; col++) {
      if (col == skipCol)
        continue;

      output[outputRow, outputCol] = input[row, col];

      outputCol++;
    }
    outputRow++;
  }

  return output;
}

矩阵逆矩阵

矩阵逆函数消除了一个函数，并通过检查矩阵秩来确定所有列是否唯一。另一个函数通过检查行列式来检查矩阵是否可逆。

/// <summary>
/// Returns a value indicates whether the specified matrix is invertible. 
/// Internally uses array determinant.
/// </summary>
public static bool IsInvertible(double[,] input) {
  var rowCount = input.GetLength(0);
  var colCount = input.GetLength(1);

  if (rowCount != colCount)
    return false;

  return Determinant(input) != 0;
}

/// <summary>
/// Calculates the inverse of a matrix. Returns null if non-invertible.
/// </summary>
public static double[,] Invert(double[,] matrix) {
  var rowCount = matrix.GetLength(0);
  var colCount = matrix.GetLength(1);

  if (rowCount != colCount)
    throw new InvalidOperationException(Properties.Resources.Exception_RequiredSquareMatrix);

  var newMatrix = ConcatHorizontally(matrix, CreateIdentityMatrix(rowCount));

  var result = Eliminate(newMatrix, MatrixReductionForm.ReducedRowEchelonForm, rowCount);
  if (result.Rank != colCount)
    return null;

  return result.Solution;
}

矩阵乘法

有四种可用的 Multiply 函数变体：标量乘法、矩阵乘法、不安全乘法和幂。

/// <summary>
/// Multiplies/Scales a matrix by a scalar input.
/// </summary>
public static double[,] Multiply(double scalar, double[,] input) {
  var rowCount = input.GetLength(0);
  var colCount = input.GetLength(1);

  // creating the final product matrix
  double[,] output = new double[rowCount, colCount];

  for (int row = 0; row < rowCount; row++) {
    for (int col = 0; col < colCount; col++) {
      output[row, col] = scalar * input[row, col];
    }
  }

  return output;
}

/// <summary>
/// Multiplies two matrices together.
/// </summary>
public static double[,] Multiply(double[,] matrix1, double[,] matrix2) {
  // caching matrix lengths for better performance
  var matrix1Rows = matrix1.GetLength(0);
  var matrix1Cols = matrix1.GetLength(1);
  var matrix2Rows = matrix2.GetLength(0);
  var matrix2Cols = matrix2.GetLength(1);

  // checking if product is defined
  if (matrix1Cols != matrix2Rows)
    throw new InvalidOperationException(Properties.Resources.Exception_UndefinedProduct);

  // creating the final product matrix
  double[,] output = new double[matrix1Rows, matrix2Cols];

  // looping through matrix 1 rows
  for (int matrix1_row = 0; matrix1_row < matrix1Rows; matrix1_row++) {
    // for each matrix 1 row, loop through matrix 2 columns
    for (int matrix2_col = 0; matrix2_col < matrix2Cols; matrix2_col++) {
      // loop through matrix 1 columns to calculate the dot product
      for (int matrix1_col = 0; matrix1_col < matrix1Cols; matrix1_col++) {
        output[matrix1_row, matrix2_col] += matrix1[matrix1_row, matrix1_col] * 
                                                matrix2[matrix1_col, matrix2_col];
      }
    }
  }

  return output;
}

/// <summary>
/// Multiplies two matrices together. Uses unsafe code. Better with extremely large matrices.
/// </summary>
public static double[,] MultiplyUnsafe(double[,] matrix1, double[,] matrix2) {
  // caching matrix lengths for better performance
  var matrix1Rows = matrix1.GetLength(0);
  var matrix1Cols = matrix1.GetLength(1);
  var matrix2Rows = matrix2.GetLength(0);
  var matrix2Cols = matrix2.GetLength(1);

  // checking if product is defined
  if (matrix1Cols != matrix2Rows)
    throw new InvalidOperationException(Properties.Resources.Exception_UndefinedProduct);

  // creating the final product matrix
  double[,] output = new double[matrix1Rows, matrix2Cols];

  unsafe
  {
    // fixing pointers to matrices
    fixed (
      double* pProduct = output,
      pMatrix1 = matrix1,
      pMatrix2 = matrix2) {

      int i = 0;
      // looping through matrix 1 rows
      for (int matrix1_row = 0; matrix1_row < matrix1Rows; matrix1_row++) {
        // for each matrix 1 row, loop through matrix 2 columns
        for (int matrix2_col = 0; matrix2_col < matrix2Cols; matrix2_col++) {
          // loop through matrix 1 columns to calculate the dot product
          for (int matrix1_col = 0; matrix1_col < matrix1Cols; matrix1_col++) {

            var val1 = *(pMatrix1 + (matrix1Rows * matrix1_row) + matrix1_col);
            var val2 = *(pMatrix2 + (matrix2Cols * matrix1_col) + matrix2_col);

            *(pProduct + i) += val1 * val2;
          }

          i++;
        }
      }
    }
  }

  return output;
}

/// <summary>
/// Raises an input matrix to the nth power.
/// </summary>
public static double[,] Power(double[,] input, int power) {
  if (input.GetLength(0) != input.GetLength(1))
    throw new InvalidOperationException(Properties.Resources.Exception_RequiredSquareMatrix);
  if (power < 0)
    throw new ArgumentException(nameof(power));
  if (power == 0)
    return CreateIdentityMatrix(input.GetLength(0));

  double[,] output = CreateCopy(input);
  for (int i = 0; i < power; i++) {
    output = Multiply(output, input);
  }

  return output;
}

子空间投影矩阵

以下函数从子空间创建投影矩阵。

/// <summary>
/// Creates projection matrix for the specified subspace.
/// </summary>
public static double[,] CreateProjectionMatrix(double[,] subspace) {
  var subspaceTranspose = MatrixFunctions.Transpose(subspace);
  double[,] value = MatrixFunctions.Multiply(subspaceTranspose, subspace);
  value = MatrixFunctions.Invert(value);
  value = MatrixFunctions.Multiply(value, subspaceTranspose);
  value = MatrixFunctions.Multiply(subspace, value);
  return value;
}

反射矩阵

以下函数创建 2D 和 3D 反射矩阵

public static double[,] Create2DReflectionMatrix(MatrixAxis axis) {
  if (axis == MatrixAxis.Z)
    throw new InvalidOperationException(Properties.Resources.Exception_InvalidAxis);

  var output = CreateIdentityMatrix(2);

  if (axis == MatrixAxis.X)
    output[1, 1] *= -1;
  else  if (axis == MatrixAxis.Y )
    output[0, 0] *= -1;

  return output;
}

public static double[,] Create3DReflectionMatrix(Matrix3DReflectionPlane axis) {
  var output = CreateIdentityMatrix(4);

  if (axis == Matrix3DReflectionPlane.XY)
    output[2, 2] *= -1;

  else if (axis == Matrix3DReflectionPlane.YZ)
    output[0,0] *= -1;

  else if (axis == Matrix3DReflectionPlane.ZX)
    output[1,1] *= -1;

  return output;
}

旋转矩阵

以下函数创建 2D 和 3D 旋转矩阵。它们接受角度、单位（弧度/度）和方向（顺时针/逆时针）

/// <summary>
/// Creates 2-dimensional rotation matrix using the specified angle.
/// </summary>
public static double[,] Create2DRotationMatrix
       (double angle, AngleUnit unit, MatrixRotationDirection direction) {
  // sin and cos accept only radians

  double angleRadians = angle;
  if (unit == AngleUnit.Degrees)
    angleRadians = Converters.DegreesToRadians(angleRadians);

  double[,] output = new double[2, 2];

  output[0, 0] = Math.Cos(angleRadians);
  output[1, 0] = Math.Sin(angleRadians);
  output[0, 1] = Math.Sin(angleRadians);
  output[1, 1] = Math.Cos(angleRadians);

  if (direction == MatrixRotationDirection.Clockwise)
    output[1, 0] *= -1;
  else
    output[0, 1] *= -1;

  return output;
}

/// <summary>
/// Creates 2-dimensional rotation matrix using the specified angle and axis.
/// </summary>
public static double[,] Create3DRotationMatrix
     (double angle, AngleUnit unit, MatrixAxis axis) {
  // sin and cos accept only radians

  double angleRadians = angle;
  if (unit == AngleUnit.Degrees)
    angleRadians = Converters.DegreesToRadians(angleRadians);

  double[,] output = new double[3, 3];

  if (axis == MatrixAxis.X) {
    output[0, 0] = 1;
    output[1, 1] = Math.Cos(angleRadians);
    output[2, 1] = Math.Sin(angleRadians);
    output[1, 2] = -1 * Math.Sin(angleRadians);
    output[2, 2] = Math.Cos(angleRadians);
  } else if (axis == MatrixAxis.Y) {
    output[1, 1] = 1;
    output[0, 0] = Math.Cos(angleRadians);
    output[2, 0] = -1 * Math.Sin(angleRadians);
    output[0, 2] = Math.Sin(angleRadians);
    output[2, 2] = Math.Cos(angleRadians);
  } else if (axis == MatrixAxis.Z) {
    output[2, 2] = 1;
    output[0, 0] = Math.Cos(angleRadians);
    output[1, 0] = Math.Sin(angleRadians);
    output[0, 1] = -1 * Math.Sin(angleRadians);
    output[1, 1] = Math.Cos(angleRadians);
  }

  return output;
}

错切矩阵

以下函数使用给定轴上的特定因子创建错切矩阵。

public static double[,] Create2DShearingMatrix(double factor, MatrixAxis axis) {
  if (axis == MatrixAxis.Z)
    throw new InvalidOperationException(Properties.Resources.Exception_InvalidAxis);

  var output = CreateIdentityMatrix(2);

  if (axis ==  MatrixAxis.X)
    output[0, 1] = factor;

  else if (axis == MatrixAxis.Y )
    output[1, 0] = factor;

  return output;
}

public static double[,] Create3DShearingMatrix(double factor, MatrixAxis axis) {
  var output = CreateIdentityMatrix(4);

  if (axis == MatrixAxis.X) {
    output[1, 0] = factor;
    output[2, 0] = factor;
  } else if (axis == MatrixAxis.Y) {
    output[0, 1] = factor;
    output[2, 1] = factor;
  } else if (axis == MatrixAxis.Z) {
    output[0, 2] = factor;
    output[1, 2] = factor;
  }

  return output;
}

向量

以下列出了此库中涵盖的一些主要向量函数。我们正在查看原始函数。包装器稍后介绍。

夹角

此函数计算两个向量之间的夹角。它依赖于另外两个函数，DotProduct 和 GetMagnitude。

public static double AngleBetween(double[] vector1, double[] vector2, AngleUnit unit) {
  if (vector1.Length != vector2.Length)
    throw new InvalidOperationException(Properties.Resources.Exception_DimensionsMismatch);

  var dotProduct = DotProduct(vector1, vector2);
  var lengthProduct = GetMagnitude(vector1) * GetMagnitude(vector2);

  var result = Math.Acos(dotProduct / lengthProduct);
  if (unit == AngleUnit.Degrees)
    result = Converters.RadiansToDegrees(result);

  return result;
}

规范化

以下代码演示了两个操作：归一化（即，转换为单位向量）和模长（长度）计算器。

/// <summary>
/// Normalizes a vector.
/// </summary>
public static double[] ToUnitVector(double[] input) {
  var length = input.Length;

  double[] output = new double[length];
  var coeffecient = 1.0 / GetMagnitude(input);

  for (int i = 0; i < length; i++) {
    output[i] = coeffecient * input[i];
  }

  return output;
}

public static double GetMagnitude(double[] input) {
  double val = 0;

  for (int i = 0; i < input.Length; i++)
    val += input[i] * input[i];

  val = Math.Sqrt(val);

  return val;
}

点积和叉积

以下代码演示了缩放、点积和叉积。

public static double DotProduct(double[] vector1, double[] vector2) {
  if (vector1.Length != vector2.Length)
    throw new InvalidOperationException(Properties.Resources.Exception_DimensionsMismatch);

  double product = 0;
  for (int i = 0; i < vector1.Length; i++)
    product += vector1[i] * vector2[i];

  return product;
}

public static double[] Scale(double scalar, double[] vector) {
  double[] product = new double[vector.Length];

  for (int i = 0; i < vector.Length; i++) {
    product[i] = scalar * vector[i];
  }

  return product;
}

public static double[] CrossProduct(double[] vector1, double[] vector2) {
  int length = 3;

  if (length != vector1.Length || length != vector2.Length)
    throw new InvalidOperationException(Properties.Resources.Exception_3DRequired);

  double[,] matrix = new double[length, length];
  for (int row = 0; row < length; row++) {
    for (int col = 0; col < length; col++) {
      if (row == 0)
        matrix[row, col] = 1;
      else if (row == 1)
        matrix[row, col] = vector1[col];
      else if (row == 2)
        matrix[row, col] = vector2[col];
    }
  }

  return MatrixFunctions.CrossProduct(matrix);
}

Projection

使用以下函数处理向量在另一个向量上的投影

public static double ProjectionFactor(double[] vector1, double[] vector2) {
  return DotProduct(vector1, vector2) / DotProduct(vector2, vector2);
}
public static double[] Projection(double[] vector1, double[] vector2) {
  var factor = ProjectionFactor(vector1, vector2);
  return Scale(factor, vector2);
}

包装器

使用原始 Array 对象（尤其是多维数组）相当麻烦。为了使其更容易，我们将这些原始函数包装到两个包装类 Matrix 和 Vector 中。这两个类都实现了 ICloneable 和 IEquatable。两者都有索引器和直接访问其内部表示（即，数组）的方式。

属性和方法

下图显示了两个类的导出成员。

运算符

我们为这两个类重写了一些运算符，以便能够以干净清晰的方式访问它们的操作。这是 Matrix 的重写运算符列表

public static Matrix operator +(Matrix m) { return m; }
public static Matrix operator -(Matrix m) { return m.Scale(-1); }
public static Matrix operator +(Matrix a, Matrix b) { return a.Add(b); }
public static Matrix operator -(Matrix a, Matrix b) { return a.Subtract(b); }
public static Matrix operator +(Matrix a, double[,] b) { return a.Add(b); }
public static Matrix operator -(Matrix a, double[,] b) { return a.Subtract(b); }
public static Matrix operator *(double a, Matrix m) { return m.Scale(a); }
public static Matrix operator *(Matrix a, Matrix b) { return a.Multiply(b); }
public static Matrix operator *(Matrix a, double[,] b) { return a.Multiply(b); }
public static bool operator ==(Matrix a, Matrix b) {
  if ((a as object) == null || (b as object) == null)
    return false;
  return a.Equals(b);
}
public static bool operator !=(Matrix a, Matrix b) {
  if ((a as object) == null || (b as object) == null)
    return false;
  return a.Equals(b) == false;
}
public static Matrix operator ^(Matrix a, int power) { return a.Power(power); }
public static implicit operator double[,] (Matrix m) { return m.InnerMatrix; }
public static explicit operator Matrix(double[,] m) { return new Linears.Matrix(m); }

最后的话

如前所述，该代码使用 MIT 许可证。因此，请在没有任何限制、约束或保证的情况下使用该代码，并随时转发您的意见、建议和更正。

历史

• 2020 年 3 月 12 日：初始版本