写出上述程序不用E
时间: 2023-12-18 08:02:35 浏览: 27
不使用第三方库的话,可以手写QR分解和矩阵运算的函数。下面是一个简单的示例代码,实现了QR分解和解超定方程组的最小二乘解。
```c++
#include <iostream>
#include <vector>
#include <cmath>
using namespace std;
// 定义矩阵类
class Matrix {
public:
Matrix(int row, int col) : row_(row), col_(col) {
data_.resize(row * col);
}
double& operator() (int i, int j) {
return data_[i * col_ + j];
}
double operator() (int i, int j) const {
return data_[i * col_ + j];
}
Matrix transpose() const {
Matrix T(col_, row_);
for (int i = 0; i < row_; ++i) {
for (int j = 0; j < col_; ++j) {
T(j, i) = (*this)(i, j);
}
}
return T;
}
Matrix operator* (const Matrix& other) const {
Matrix C(row_, other.col_);
for (int i = 0; i < row_; ++i) {
for (int j = 0; j < other.col_; ++j) {
double s = 0;
for (int k = 0; k < col_; ++k) {
s += (*this)(i, k) * other(k, j);
}
C(i, j) = s;
}
}
return C;
}
void qr(Matrix& Q, Matrix& R) const {
Q = Matrix(row_, col_);
R = Matrix(col_, col_);
for (int j = 0; j < col_; ++j) {
// 计算第 j 列的范数
double norm = 0;
for (int i = j; i < row_; ++i) {
norm = hypot(norm, (*this)(i, j));
}
if (norm == 0) {
// 如果范数为 0,说明第 j 列已经为零向量,不需要进行变换
R(j, j) = 0;
continue;
}
// 将第 j 列缩放为单位向量
for (int i = j; i < row_; ++i) {
(*this)(i, j) /= norm;
}
// 计算第 j 列对应的 Householder 变换
double s = 0;
for (int i = j; i < row_; ++i) {
s += (*this)(i, j) * (*this)(i, j);
}
double beta = 2 / s;
for (int k = j + 1; k < col_; ++k) {
double s = 0;
for (int i = j; i < row_; ++i) {
s += (*this)(i, j) * (*this)(i, k);
}
double t = beta * s;
for (int i = j; i < row_; ++i) {
(*this)(i, k) -= t * (*this)(i, j);
}
}
// 记录 Q 和 R 的值
for (int i = 0; i < row_; ++i) {
Q(i, j) = (*this)(i, j);
}
for (int k = j; k < col_; ++k) {
R(j, k) = 0;
for (int i = j; i < row_; ++i) {
R(j, k) += Q(i, j) * (*this)(i, k);
}
}
}
}
Matrix solve(const Matrix& b) const {
Matrix Q, R;
qr(Q, R);
Matrix Qt_b = Q.transpose() * b;
Matrix x(col_, 1);
for (int i = col_ - 1; i >= 0; --i) {
double s = 0;
for (int j = i + 1; j < col_; ++j) {
s += R(i, j) * x(j, 0);
}
x(i, 0) = (Qt_b(i, 0) - s) / R(i, i);
}
return x;
}
private:
int row_, col_;
vector<double> data_;
};
int main() {
// 读入数据,假设有 n 个数据点和 m 个未知数
int n = 5, m = 3;
Matrix A(n, m);
Matrix b(n, 1);
for (int i = 0; i < n; ++i) {
double x, y, z;
// 读入第 i 个数据点的 x, y, z 坐标
// ...
// 计算第 i 行系数矩阵和结果向量的值
A(i, 0) = x; A(i, 1) = y; A(i, 2) = 1;
b(i, 0) = z;
}
// 解方程 Ax = b 的最小二乘解
Matrix x = A.solve(b);
// 输出结果
for (int i = 0; i < m; ++i) {
cout << "x[" << i << "] = " << x(i, 0) << endl;
}
return 0;
}
```
在这个示例代码中,我们手写了矩阵类、QR分解和解超定方程组的最小二乘解函数。在读入数据后,我们调用 `A.solve(b)` 函数求解方程组的最小二乘解。最后输出结果即可。
相关推荐
![rar](https://img-home.csdnimg.cn/images/20210720083606.png)
![zip](https://img-home.csdnimg.cn/images/20210720083736.png)
![zip](https://img-home.csdnimg.cn/images/20210720083736.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)