从零实现矩阵求逆与线性方程组求解：深入理解高斯消元法

1. 引言

在机器学习和科学计算中，矩阵求逆和线性方程组求解是基础而重要的操作。虽然 NumPy 等库提供了高效的实现，但理解其底层原理对于深入掌握线性代数至关重要。本文将带你从零实现矩阵求逆和线性方程组求解算法，通过代码实践深入理解高斯消元法的原理。

2. 矩阵求逆的实现

2.1 算法原理

矩阵求逆的核心思想是高斯-约旦消元法。对于一个 n×n 的方阵 A，我们构造一个增广矩阵 [A|I]，其中 I 是 n×n 的单位矩阵。通过对增广矩阵进行行变换，将 A 部分化为单位矩阵，此时右侧的 I 部分就变成了 A 的逆矩阵。

2.2 代码实现

import numpy as np

def matrix_inverse(a):
"""计算方阵的逆矩阵"""
n = a.shape[0]
if a.shape[0] != a.shape[1]:
raise ValueError("非方阵，无法计算逆矩阵。")

# 构造增广矩阵 [A|I]
a_identity = np.hstack((a, np.eye(a.shape[0])))

# 进行高斯-约旦消元
for i in range(n):
# 列主元选择：找到第i列绝对值最大的行
max_row = i
for j in range(i, n):
if np.abs(a_identity[j, i]) > np.abs(a_identity[max_row, i]):
max_row = j

# 交换当前行与主元行
if max_row != i:
a_identity[[i, max_row], :] = a_identity[[max_row, i], :]

pivot = a_identity[i, i]
if abs(pivot) < 1e-12:
raise ValueError("矩阵奇异，不存在逆矩阵")

# 将主元行归一化
a_identity[i, :] = a_identity[i, :] / pivot

# 消去其他行的第i列元素
for j in range(n):
if j != i:
factor = a_identity[j, i]
a_identity[j, :] = a_identity[j, :] – a_identity[i, :] * factor

# 返回逆矩阵部分
return a_identity[:, n:]

# 测试代码
if __name__ == "__main__":
# 测试矩阵求逆
a = np.array([[1, 2, 3], [4, 8, 6], [7, 6, 9]])
print("原始矩阵 A:")
print(a)
print("\\nNumPy 计算的逆矩阵:")
print(np.linalg.inv(a))
print("\\n自定义实现的逆矩阵:")
custom_inv = matrix_inverse(a)
print(custom_inv)

# 验证结果
if np.allclose(np.linalg.inv(a), custom_inv, rtol=1e-10):
print("\\n✅ 自定义逆矩阵计算正确")
else:
print("\\n❌ 自定义逆矩阵计算错误")

2.3 关键优化点

3. 矩阵乘法的实现

3.1 基础矩阵乘法

def matrix_dot(a, b):
"""实现矩阵乘法，兼容向量输入"""
# 兼容一维向量输入
if len(a.shape) == 1:
a = a.reshape(1, –1)
if len(b.shape) == 1:
b = b.reshape(–1, 1)

# 检查维度匹配
ka = a.shape[1]
kb = b.shape[0]
if ka != kb:
raise ValueError(f"矩阵维度不匹配：a的列{ka}，b的行{kb}，无法相乘")

# 三重循环实现矩阵乘法
result = np.zeros((a.shape[0], b.shape[1]))
for i in range(a.shape[0]):
for j in range(b.shape[1]):
s = 0.0
for k in range(a.shape[1]):
s += a[i, k] * b[k, j]
result[i, j] = s

return result

# 测试矩阵乘法
def test_matrix_dot():
a = np.array([[1, 2], [3, 4]])
b = np.array([[5, 6], [7, 8]])

print("矩阵 A:")
print(a)
print("\\n矩阵 B:")
print(b)
print("\\nNumPy 计算结果:")
print(np.dot(a, b))
print("\\n自定义实现结果:")
print(matrix_dot(a, b))

# 验证向量乘法
v1 = np.array([1, 2, 3])
v2 = np.array([4, 5, 6])
print("\\n向量点积验证:")
print("NumPy:", np.dot(v1, v2))
print("自定义:", matrix_dot(v1, v2))

4. 线性方程组求解

4.1 基于逆矩阵的求解

def linear_equations_solve(a, b):
"""使用逆矩阵法求解线性方程组 Ax = b"""
# 记录原始b是否为一维向量
is_1d = len(b.shape) == 1

# 计算逆矩阵
inverse_a = matrix_inverse(a)

# 计算解：x = A⁻¹b
res = matrix_dot(inverse_a, b)

# 如果原始b是一维向量，将结果转换为一维
if is_1d:
res = res.flatten()

return res

# 测试线性方程组求解
def test_linear_solve():
# 测试用例1：一般矩阵
a = np.array([[1, 2, 3], [4, 8, 6], [7, 6, 9]])
b = np.array([1, 2, 3])

print("测试用例1 – 一般矩阵")
print("系数矩阵 A:")
print(a)
print("\\n常数向量 b:", b)

print("\\nNumPy 求解结果:")
np_solution = np.linalg.solve(a, b)
print(np_solution)

print("\\n自定义求解结果:")
custom_solution = linear_equations_solve(a, b)
print(custom_solution)

if np.allclose(np_solution, custom_solution, rtol=1e-10):
print("\\n✅ 线性方程组求解正确")
else:
print("\\n❌ 线性方程组求解错误")

# 测试用例2：单位矩阵
print("\\n\\n测试用例2 – 单位矩阵")
a_identity = np.eye(3)
b_vec = np.array([1, 2, 3])

print("系数矩阵 A (单位矩阵):")
print(a_identity)
print("\\n常数向量 b:", b_vec)

print("\\nNumPy 求解结果:")
print(np.linalg.solve(a_identity, b_vec))

print("\\n自定义求解结果:")
print(linear_equations_solve(a_identity, b_vec))

5. 算法复杂度分析

5.1 时间复杂度

5.2 空间复杂度

• 矩阵求逆：O(n²)，需要存储增广矩阵
• 矩阵乘法：O(n²)，需要存储结果矩阵
• 整体算法：O(n²)

6. 数值稳定性考虑

6.1 主元选择的重要性

def test_pivot_importance():
"""演示主元选择对数值稳定性的影响"""
# 创建一个条件数较大的矩阵
a = np.array([[1e-10, 1], [1, 1]])
b = np.array([1, 2])

print("条件数较大的矩阵:")
print("A =", a)
print("b =", b)

try:
# 不使用主元选择（简化版本）
print("\\n不使用主元选择的结果:")
# 这里可以对比不同实现
except Exception as e:
print(f"计算失败: {e}")

print("\\n使用列主元选择的结果:")
solution = linear_equations_solve(a, b)
print("x =", solution)

6.2 浮点数精度处理

在实现中需要注意：

7. 性能优化建议

7.1 向量化优化

def matrix_dot_optimized(a, b):
"""优化版的矩阵乘法"""
if len(a.shape) == 1:
a = a.reshape(1, –1)
if len(b.shape) == 1:
b = b.reshape(–1, 1)

# 使用 @ 运算符或 np.dot 的优化版本
return a @ b # 或者 np.dot(a, b)

def matrix_inverse_optimized(a):
"""优化版的矩阵求逆"""
n = a.shape[0]
a_identity = np.hstack((a, np.eye(n)))

for i in range(n):
# 列主元选择（向量化）
max_row = i + np.argmax(np.abs(a_identity[i:, i]))
if max_row != i:
a_identity[[i, max_row]] = a_identity[[max_row, i]]

pivot = a_identity[i, i]
if abs(pivot) < 1e-12:
raise ValueError("矩阵奇异")

# 归一化当前行
a_identity[i] /= pivot

# 消去其他行（向量化）
for j in range(n):
if j != i:
a_identity[j] -= a_identity[i] * a_identity[j, i]

return a_identity[:, n:]

7.2 使用 BLAS/LAPACK

对于生产环境，建议：

8. 总结

本文从零实现了矩阵求逆和线性方程组求解的核心算法，重点包括：

虽然这些实现不如 NumPy 优化，但理解其原理对于深入学习线性代数和数值计算至关重要。在实际应用中，建议使用经过高度优化的数值计算库，但在学习和教学场景中，手动实现这些算法是极好的练习。

9. 扩展练习

# 扩展：验证逆矩阵性质
def verify_inverse_properties(a):
"""验证逆矩阵的基本性质"""
inv_a = matrix_inverse(a)

# 性质1: A × A⁻¹ = I
product1 = matrix_dot(a, inv_a)
identity1 = np.eye(a.shape[0])

# 性质2: A⁻¹ × A = I
product2 = matrix_dot(inv_a, a)
identity2 = np.eye(a.shape[0])

print(f"A × A⁻¹ 接近单位矩阵: {np.allclose(product1, identity1)}")
print(f"A⁻¹ × A 接近单位矩阵: {np.allclose(product2, identity2)}")

return inv_a

通过本文的学习，你应该对矩阵求逆和线性方程组求解有了更深入的理解。继续探索数值线性代数的世界，你会发现更多有趣的算法和应用！

附

以上为AI基于以下测试验证代码生成，可以作为参考。

import numpy as np

def matrix_inverse(a):
n = a.shape[0]
if a.shape[0] != a.shape[1]:
raise ValueError("非方阵，无法计算逆矩阵。")

a_identity = np.hstack((a, np.eye(a.shape[0])))
#先把扩展为n*2n的矩阵
# a_identity = np.zeros((n, 2*n), dtype=float)
# a_identity[:, :n] = a
# #先拼接基础矩阵
# for i in range(n):
# for j in range(n, a_identity.shape[1]):
# if i == j – n:
# a_identity[i, j] = 1
# else:
# a_identity[i, j] = 0
# 进行高斯消元法
for i in range(n):
#找到本列绝对值最大的行
max_row = i
for j in range(i, n): # 从i行向下开始查找绝对值最大的行，上面的是已经交换过的
if np.abs(a_identity[j, i]) > np.abs(a_identity[max_row, i]):
max_row = j
temp = a_identity[i, :].copy()
a_identity[i, :] = a_identity[max_row, :]
a_identity[max_row, :] = temp
# a_identity[[i, max_row], :] = i_identity[[max_row, i], :] # 替换，交换i行和max_row行

pivot = a_identity[i, i]
if abs(pivot) < 1e-12:
raise ValueError("矩阵奇异，不存在逆矩阵")

# 对i行进行消元，将i行的第i列元素设为1
for j in range(a_identity.shape[1]):
a_identity[i, j] = a_identity[i, j] / pivot
# a_identity[i, :] = a_identity[i,:]/a_identity[i, i] # 替换，对i行进行消元，将i行的第i列元素设为1

# 对其他行进行消元，将其他行的第i列元素设为0
for j in range(n):
if j != i:
factor = a_identity[j, i] # ✅提前取出系数
for k in range(a_identity.shape[1]):
a_identity[j, k] = a_identity[j, k] – a_identity[i, k] * factor
# a_identity[j, :] = a_identity[j,:] – a_identity[i,:] * a_identity[j, i] # 替换 对其他行进行消元，将其他行的第i列元素设为0
return a_identity[:, n:]

def matrix_dot(a, b):
# 兼容一维向量输入
if len(a.shape) == 1:
a = a.reshape(1, –1)
if len(b.shape) == 1:
b = b.reshape(–1, 1)

ka = a.shape[1]
kb = b.shape[0]

if ka != kb:
raise ValueError(f"矩阵维度不匹配：a的列{ka}，b的行{kb}，无法相乘")

result = np.zeros((a.shape[0], b.shape[1]))
for i in range(a.shape[0]):
for j in range(b.shape[1]):
s = 0.0
for k in range(a.shape[1]):
s += a[i, k] * b[k,j]
result[i,j] = s
return result

def linear_quations_solve(a, b):
# 求解线性方程组
# 先求a的逆矩阵
# 再用逆矩阵乘以b
# 即：x = a^{-1} * b
# 记录原始b是否一维
is_1d = len(b.shape) == 1
inverse_a = matrix_inverse(a)
res = matrix_dot(inverse_a, b)
if is_1d:
res = res.flatten() # 将结果转换为一维向量，ravel()方法会改变原始数组的形状
return res

if __name__ == "__main__":
a = np.array([[1, 2, 3], [4, 8, 6], [7, 6, 9]])
print(a)
print("a的逆矩阵为：")
print(np.linalg.inv(a))
print("自定义逆矩阵为：")
print(matrix_inverse(a))
if np.allclose(np.linalg.inv(a), matrix_inverse(a)): # 不能直接用==，因为浮点数有精度问题
print("自定义逆矩阵计算正确")
else:
print("自定义逆矩阵计算错误")

#求解线性方程组
b = np.array([1, 2, 3])
x = np.linalg.solve(a, b)
print("线性方程组的解为：")
print(x)
print("自定义线性方程组的解为：")
print(linear_quations_solve(a, b))

if np.allclose(np.linalg.solve(a, b), linear_quations_solve(a, b)):
print("自定义线性方程组求解正确")
else:
print("自定义线性方程组求解错误")

#求解线性方程组
a1 = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
print(np.linalg.inv(a1))
b1 = np.array([1, 2, 3])
x1 = np.linalg.solve(a1, b1)
print("线性方程组的解为：")
print(x1)

print("矩阵的逆矩阵为：")
print(matrix_inverse(a1))