C#，码海拾贝（40）——求解“线性最小二乘问题”的“豪斯荷尔德Householder变换法”之C#源代码

// C#代码示例using System;

class Householder{
static void Main()
{
// 定义矩阵A和向量b double[,] A = {{1,2,3}, {4,5,6}, {7,8,9}};
double[] b = {10,11,12};

// 调用Householder变换函数 double[] x = HouseholderSolve(A, b);

// 输出结果 Console.WriteLine("解向量x为：");
for (int i =0; i < x.Length; i++)
{
Console.WriteLine(x[i]);
}
}

// Householder变换函数 static double[] HouseholderSolve(double[,] A, double[] b)
{
int n = A.GetLength(0);
double[] x = new double[n];

for (int k =0; k < n -1; k++)
{
double[] v = new double[n - k];
double[] u = new double[n - k];

double sigma =0;
for (int i = k; i < n; i++)
{
sigma += A[i, k] * A[i, k];
}

double a = Math.Sqrt(sigma);
if (A[k, k] >0)
{
a = -a;
}

v[0] = A[k, k] - a;
for (int i = k +1; i < n; i++)
{
v[i - k] = A[i, k];
}

double beta =2 / (v[0] * v[0] + sigma);
for (int i =0; i < n - k; i++)
{
u[i] = v[i] * beta;
}

for (int j = k; j < n; j++)
{
double gamma =0;
for (int i =0; i < n - k; i++)
{
gamma += u[i] * A[j, k + i];
}
for (int i =0; i < n - k; i++)
{
A[j, k + i] -= gamma * u[i];
}
}

double[] y = new double[n - k];
for (int i =0; i < n - k; i++)
{
double temp =0;
for (int j =0; j < n - k; j++)
{
temp += A[k + i, k + j] * b[k + j];
}
y[i] = temp;
}

for (int i =0; i < n - k; i++)
{
double temp =0;
for (int j =0; j < n - k; j++)
{
temp += A[k + i, k + j] * y[j];
}
b[k + i] -= temp;
}
}

x[n -1] = b[n -1] / A[n -1, n -1];
for (int i = n -2; i >=0; i--)
{
double temp =0;
for (int j = i +1; j < n; j++)
{
temp += A[i, j] * x[j];
}
x[i] = (b[i] - temp) / A[i, i];
}

return x;
}
}