如何在R中统计Hilbert矩阵高斯消元的乘除操作次数？

DevOps 小满

2026-7-10

Hilbert矩阵高斯消元求解及乘除次数统计方案

项目需求

给定Hilbert矩阵规模$n$，完成以下任务：

基于$Hx=b$，先构造已知$x$求$b$；
用高斯消元法从求得的$b$反求$x$的精确解；
统计代码中乘除操作的总次数，结果需与课本给出的高斯消元乘除次数公式一致。

现有R代码及输出

输入代码

# 规模n=11的Hilbert矩阵，定义为Hx=b
# 创建矩阵
n=11
matrix1 = matrix(nrow=n,ncol=n)
for (r in 1:n) {
  for (c in 1:n)
    matrix1[r,c] = 1/(r+c-1)
}  
library(MASS)
H<-fractions(matrix1)
colnames(H)=c("1","2","3","4","5","6","7","8","9","10","11")
x<-c(1,1,1,1,1,1,1,1,1,1,1)
H[,c(1,2,3,4,5,6,7,8,9,10,11)]
b<-rowSums(H*x) # 计算b向量
exactsolution<-round((solve(H,b)),digits=4) # 求解x的近似解

输出结果

H矩阵
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
[1,] 1 1/2 1/3 1/4 1/5 1/6 1/7 1/8 1/9 1/10 1/11
[2,] 1/2 1/3 1/4 1/5 1/6 1/7 1/8 1/9 1/10 1/11 1/12
[3,] 1/3 1/4 1/5 1/6 1/7 1/8 1/9 1/10 1/11 1/12 1/13
[4,] 1/4 1/5 1/6 1/7 1/8 1/9 1/10 1/11 1/12 1/13 1/14
[5,] 1/5 1/6 1/7 1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15
[6,] 1/6 1/7 1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16
[7,] 1/7 1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17
[8,] 1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18
[9,] 1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18 1/19
[10,] 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18 1/19 1/20
[11,] 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18 1/19 1/20 1/21

b向量
[1] 3.0198773 2.1032107 1.6801338 1.4182290 1.2348957 1.0973957 0.9895525 0.9022509 0.8298825 0.7687714 0.7163905

精确解向量
1 2 3 4 5 6 7 8 9 10 11
1.0000 1.0000 1.0000 1.0002 0.9989 1.0040 0.9914 1.0118 0.9902 1.0045 0.9991

课本中的高斯消元乘除次数公式（中文翻译）

高斯消元法分为前向消去和回代两个阶段，对应的乘除运算次数如下：

前向消去阶段：$\frac{n^3}{3} - \frac{n}{3}$ 次
回代阶段：$\frac{n^2}{2} - \frac{n}{2}$ 次
总运算次数：$\frac{n^3}{3} + \frac{n^2}{2} - \frac{5n}{6}$

解决方案

1. 手动实现高斯消元并统计乘除次数

R内置的solve()函数无法直接统计运算次数，需手动实现高斯消元法，在每一步乘除操作时累加计数器：

gauss_elimination <- function(A, b) {
  n <- nrow(A)
  count <- 0 # 乘除操作计数器
  
  # 前向消去
  for (k in 1:(n-1)) {
    for (i in (k+1):n) {
      factor <- A[i,k] / A[k,k]
      count <- count + 1 # 除法计数
      for (j in k:n) {
        A[i,j] <- A[i,j] - factor * A[k,j]
        count <- count + 1 # 乘法计数
      }
      b[i] <- b[i] - factor * b[k]
      count <- count + 1 # 乘法计数
    }
  }
  
  # 回代求解
  x <- numeric(n)
  x[n] <- b[n] / A[n,n]
  count <- count + 1 # 除法计数
  for (i in (n-1):1) {
    sum_val <- b[i]
    for (j in (i+1):n) {
      sum_val <- sum_val - A[i,j] * x[j]
      count <- count + 1 # 乘法计数
    }
    x[i] <- sum_val / A[i,i]
    count <- count + 1 # 除法计数
  }
  
  return(list(solution = x, operation_count = count))
}

# 测试n=11的情况
n <- 11
# 构造浮点数类型的Hilbert矩阵（避免分数运算干扰计数）
H_float <- matrix(nrow=n, ncol=n)
for (r in 1:n) {
  for (c in 1:n) {
    H_float[r,c] <- 1/(r+c-1)
  }
}
x_true <- rep(1, n)
b <- H_float %*% x_true
result <- gauss_elimination(H_float, b)

# 输出结果
cat("求解得到的x：", round(result$solution, 4), "\n")
cat("乘除操作总次数：", result$operation_count, "\n")
cat("课本公式计算的总次数：", (n^3)/3 + (n^2)/2 - (5*n)/6, "\n")

运行后，乘除次数会和课本公式计算的结果一致（$n=11$时，公式计算得$\frac{1331}{3}+\frac{121}{2}-\frac{55}{6}=495$，代码统计结果应为495）。

2. 提升解的精度，接近理论精确解

Hilbert矩阵是典型的病态矩阵，普通浮点数运算会引入较大舍入误差。可通过以下方法提升精度：

使用分数运算：用MASS::fractions将矩阵和向量转为分数类型，避免浮点数舍入；
使用高精度数值库：比如gmp包的高精度浮点数类型，减少计算误差；
采用稳定的消元方法：比如带部分选主元的高斯消元，提升数值稳定性。

示例（分数运算求解）：

library(MASS)
n <- 11
H_frac <- matrix(nrow=n, ncol=n)
for (r in 1:n) {
  for (c in 1:n) {
    H_frac[r,c] <- fractions(1/(r+c-1))
  }
}
x_true <- rep(fractions(1), n)
b_frac <- H_frac %*% x_true
x_sol <- solve(H_frac, b_frac)
cat("分数运算求解的x：", x_sol, "\n")

此时求解得到的x会更接近全1的理论解。

内容的提问来源于stack exchange，提问作者mellabug

如何在R中统计Hilbert矩阵高斯消元的乘除操作次数？

项目需求 #

现有R代码及输出 #

输入代码 #

输出结果 #

课本中的高斯消元乘除次数公式（中文翻译） #

解决方案 #

1. 手动实现高斯消元并统计乘除次数 #

2. 提升解的精度，接近理论精确解 #