PLU 分解的矩阵逆
Matrix inverse from PLU decomposition
我正在尝试从 plu 分解计算逆矩阵,但我得到的结果是错误的。使用调试器,我找不到任何可以责备的步骤,所以这里是反转函数的代码:
template<class T>
Matrix<T> Matrix<T>::inv() const {
std::tuple<Matrix<T>, Matrix<T>, Matrix<T>> PLU(this->decomp_PLU());
Matrix<T> P(std::get<0>(PLU));
Matrix<T> L(std::get<1>(PLU));
Matrix<T> U(std::get<2>(PLU));
if (!this->lines() == this->cols()) { throw std::logic_error("Matrix must be square"); }
unsigned N = this->lines();
Matrix<T> Y(N, N);
for (unsigned i = 0; i < N; i++) {
Matrix<T> B(Matrix<T>::gen_col(P.transp().data[i])); // B is the ith column vector of P
Y[i] = Matrix<T>::solve_climb(L, B).transp().data[0]; // Compute LY=P for each column vector
Matrix<T> conf(L.dot(Matrix<T>::gen_col(Y[i])));
if (!conf.allclose(B, 1e-9, 1e-9)) {
throw std::logic_error("WHYY");
}
}
Y = Y.transp();
Matrix<T> X(N, N);
for (unsigned i = 0; i < N; i++) {
Matrix<T> B(Matrix<T>::gen_col(Y.transp().data[i])); // B is the ith column vector of Y
X[i] = Matrix<T>::solve_descent(U, B).transp().data[0]; // Compute UX=Y for each column vector
Matrix<T> conf(U.dot(Matrix<T>::gen_col(X[i])));
if (!conf.allclose(B, 1e-9, 1e-9)) {
throw std::logic_error("WHYY");
}
}
X = X.transp();
return X;
}
函数 Matrix<T>::gen_col
从其输入生成列向量,solve_climb
/ solve_descent
分别计算求解 AX=B 的列向量,其中 A 是三角矩阵(劣或上级视情况而定)
仅供参考,当尝试检查每个向量的计算是否正确时,代码抛出逻辑错误 'WHYY'。
代码可能哪里出错的任何线索?
谢谢
编辑:完整代码是here (matrix_def.h) and here (matrix.h)
由于L是三角劣等,你应该使用solve_descent
,而不是solve_climb
。
U(三角形上级)也是一样,需要使用solve_climb
.
我正在尝试从 plu 分解计算逆矩阵,但我得到的结果是错误的。使用调试器,我找不到任何可以责备的步骤,所以这里是反转函数的代码:
template<class T>
Matrix<T> Matrix<T>::inv() const {
std::tuple<Matrix<T>, Matrix<T>, Matrix<T>> PLU(this->decomp_PLU());
Matrix<T> P(std::get<0>(PLU));
Matrix<T> L(std::get<1>(PLU));
Matrix<T> U(std::get<2>(PLU));
if (!this->lines() == this->cols()) { throw std::logic_error("Matrix must be square"); }
unsigned N = this->lines();
Matrix<T> Y(N, N);
for (unsigned i = 0; i < N; i++) {
Matrix<T> B(Matrix<T>::gen_col(P.transp().data[i])); // B is the ith column vector of P
Y[i] = Matrix<T>::solve_climb(L, B).transp().data[0]; // Compute LY=P for each column vector
Matrix<T> conf(L.dot(Matrix<T>::gen_col(Y[i])));
if (!conf.allclose(B, 1e-9, 1e-9)) {
throw std::logic_error("WHYY");
}
}
Y = Y.transp();
Matrix<T> X(N, N);
for (unsigned i = 0; i < N; i++) {
Matrix<T> B(Matrix<T>::gen_col(Y.transp().data[i])); // B is the ith column vector of Y
X[i] = Matrix<T>::solve_descent(U, B).transp().data[0]; // Compute UX=Y for each column vector
Matrix<T> conf(U.dot(Matrix<T>::gen_col(X[i])));
if (!conf.allclose(B, 1e-9, 1e-9)) {
throw std::logic_error("WHYY");
}
}
X = X.transp();
return X;
}
函数 Matrix<T>::gen_col
从其输入生成列向量,solve_climb
/ solve_descent
分别计算求解 AX=B 的列向量,其中 A 是三角矩阵(劣或上级视情况而定)
仅供参考,当尝试检查每个向量的计算是否正确时,代码抛出逻辑错误 'WHYY'。
代码可能哪里出错的任何线索?
谢谢
编辑:完整代码是here (matrix_def.h) and here (matrix.h)
由于L是三角劣等,你应该使用solve_descent
,而不是solve_climb
。
U(三角形上级)也是一样,需要使用solve_climb
.