求解线性方程组AX=B

H_yN
183 阅读4分钟

求逆矩阵AX=B

#include <eigen/Eigen/Dense>
void test01() 
{
    Eigen::Matrix3d A;
    Eigen::Vector3d b;
    A << 4, 7, 2, 
        3, 5, 1, 
        6, 8, 9;
    b << 1, 2, 3;
    if (A.determinant() != 0)   //判定是否可逆
    {
        //求逆矩阵
        Eigen::Matrix3d A_inv = A.inverse();
        //求x
        Eigen::Vector3d x = A_inv * b;
        //输出
        std::cout << "求解的x:" << x << std::endl;
    }
    else
    {
        std::cout << "不可逆,无解!" << std::endl;
    }
}

梯度下降算法求矩阵AX=B，效果很差

void test02() 
{
    // 创建一个示例线性方程组 A * X = B
    Eigen::MatrixXd A(3, 3);
    A << 4, 7, 2, 3, 5, 1, 6, 8, 9;
    Eigen::VectorXd B(3);
    B << 1, 2, 3;
    // 初始猜测值
    Eigen::VectorXd X(3);
    X << 1, 1, 1;
    // 学习率
    double learning_rate = 0.007;
    // 计算目标函数 0.5 * ||A * X - B||^2
    auto objective = [&](const Eigen::MatrixXd& A, const Eigen::VectorXd& X,
                         const Eigen::VectorXd& B) {
        Eigen::VectorXd residual = A * X - B;
        return 0.5 * residual.squaredNorm();
    };
    // 计算目标函数的梯度
    auto gradient = [&](const Eigen::MatrixXd& A, const Eigen::VectorXd& X,
                        const Eigen::VectorXd& B) {
        return A.transpose() * (A * X - B);
    };
    // 最大迭代次数和容忍误差
    int maxIter = 10000;
    double tol = 1e-6;
    int iter = 0;
    // 梯度下降迭代
    while (iter < maxIter)
    {
        // 计算当前的梯度
        Eigen::VectorXd grad = gradient(A, X, B);
        // 判断是否收敛
        if (grad.norm() < tol)
        {
            std::cout << "Converged at iteration " << iter << std::endl;
            break;
        }
        // 更新 X
        X -= learning_rate * grad;
        // 打印每隔一定迭代次数的结果
        if (iter % 100 == 0)
        {
            std::cout << "Iteration " << iter
                      << ", Objective Value: " << objective(A, X, B)
                      << std::endl;
        }
        iter++;
    }
    // 输出结果
    std::cout << "Solution X: " << X.transpose() << std::endl;

}

牛顿法求解AX=B

void test03() 
{
    Eigen::MatrixXd A(3, 3);
    Eigen::Vector3d B(3);
    A << 4, 7, 2, 3, 5, 1, 6, 8, 9;
    B << 1, 2, 3;
    Eigen::Vector3d x = Eigen::Vector3d::Ones();
    int max_iterations = 100;
    double tolerance = 1e-6;
    // 牛顿法迭代
    for (int iteration = 0; iteration < max_iterations; ++iteration)
    {
        // 计算 A * x
        Eigen::Vector3d Ax = A * x;
        // 计算误差 A * x - B
        Eigen::Vector3d error = Ax - B;
        // 计算梯度： A^T * (A * x - B)
        Eigen::Vector3d gradient = A.transpose() * error;
        // 计算海森矩阵： A^T * A
        Eigen::Matrix3d hessian = A.transpose() * A;
        // 使用海森矩阵的逆更新 x
        Eigen::Vector3d x_new = x - hessian.inverse() * gradient;
        // 判断收敛条件：如果变化小于容忍度，则停止迭代
        if ((x_new - x).norm() < tolerance)
        {
            std::cout << "Converged at iteration " << iteration + 1 << std::endl;
            break;
        }
        // 更新 x
        x = x_new;
    }
    std::cout << "Final solution x: " << x << std::endl;
}

拟牛顿法BFGS

void test04() 
{
    // 创建一个示例线性方程组 A * X = B
    Eigen::MatrixXd A(3, 3);
    A << 4, 7, 2, 3, 5, 1, 6, 8, 9;
    Eigen::VectorXd B(3);
    B << 1, 2, 3;
    // 初始猜测值
    Eigen::VectorXd X(3);
    X << 0, 0, 0;
    // 计算目标函数 0.5 * ||A * X - B||^2
    auto objective = [&](const Eigen::MatrixXd& A, const Eigen::VectorXd& X,
                         const Eigen::VectorXd& B) {
        Eigen::VectorXd residual = A * X - B;
        return 0.5 * residual.squaredNorm();
    };
    // 计算目标函数的梯度
    auto gradient = [&](const Eigen::MatrixXd& A, const Eigen::VectorXd& X,
                        const Eigen::VectorXd& B) {
        return A.transpose() * (A * X - B);
    };
    // BFGS 方法
    int n = A.cols();                         // 变量的维度
    Eigen::VectorXd grad = gradient(A, X, B); // 初始梯度
    // 初始化 H 矩阵为单位矩阵
    Eigen::MatrixXd H = Eigen::MatrixXd::Identity(n, n);
    double tol = 1e-6; // 收敛容忍度
    int maxIter = 100; // 最大迭代次数
    int iter = 0;
    while (grad.norm() > tol && iter < maxIter)
    {
        // 计算搜索方向 p = -H * grad
        Eigen::VectorXd p = -H * grad;
        // 线搜索：固定步长为1.0，可以更复杂的线搜索方法
        double alpha = 1.0;
        Eigen::VectorXd X_new = X + alpha * p;
        // 计算新的梯度
        Eigen::VectorXd grad_new = gradient(A, X_new, B);
        // 计算差值 s = X_new - X 和 y = grad_new - grad
        Eigen::VectorXd s = X_new - X;
        Eigen::VectorXd y = grad_new - grad;
        // 更新 H 矩阵
        double rho = 1.0 / y.dot(s); // rho = 1 / (y^T * s)
        Eigen::MatrixXd I = Eigen::MatrixXd::Identity(n, n);
        Eigen::MatrixXd term1 = I - rho * s * y.transpose();
        Eigen::MatrixXd term2 = I - rho * y * s.transpose();
        // H = H + (I - rho * s * y^T) * H * (I - rho * y * s^T) + rho * s * s^T
        H = term1 * H * term2 + rho * s * s.transpose();
        // 更新 X 和 grad
        X = X_new;
        grad = grad_new;
        iter++;
    }
    // 输出结果
    std::cout << "解 X: " << X.transpose() << std::endl;
}

拟牛顿法L-BFGS

void test05() 
{
    // 创建一个示例线性方程组 A * X = B
    Eigen::MatrixXd A(3, 3);
    A << 4, 7, 2, 3, 5, 1, 6, 8, 9;
    Eigen::VectorXd B(3);
    B << 1, 2, 3;
    // 初始猜测值
    Eigen::VectorXd X(3);
    X << 0, 0, 0;
    // 计算目标函数 0.5 * ||A * X - B||^2
    auto objective = [&](const Eigen::MatrixXd& A, const Eigen::VectorXd& X,
                         const Eigen::VectorXd& B) {
        Eigen::VectorXd residual = A * X - B;
        return 0.5 * residual.squaredNorm();
    };
    // 计算目标函数的梯度
    auto gradient = [&](const Eigen::MatrixXd& A, const Eigen::VectorXd& X,
                        const Eigen::VectorXd& B) {
        return A.transpose() * (A * X - B);
    };
    // L-BFGS 方法
    int n = A.cols();  // 变量的维度
    int m = 5;         // 存储的差分对的数量
    double tol = 1e-6; // 收敛容忍度
    int maxIter = 100; // 最大迭代次数
    int iter = 0;
    // 存储最近的 s 和 y
    std::vector<Eigen::VectorXd> s_list, y_list;
    // 初始梯度
    Eigen::VectorXd grad = gradient(A, X, B);
    // L-BFGS迭代
    while (grad.norm() > tol && iter < maxIter)
    {
        // 计算搜索方向
        Eigen::VectorXd q = grad;
        std::vector<double> alpha_list; // 用于存储 alpha 值
        // 计算方向 q = -H_k * grad，H_k 是通过 s 和 y 计算的近似 Hessian
        for (int i = s_list.size() - 1; i >= 0; --i)
        {
            Eigen::VectorXd s = s_list[i];
            Eigen::VectorXd y = y_list[i];
            double rho = 1.0 / y.dot(s);
            double alpha = rho * s.dot(q);
            alpha_list.push_back(alpha); // 存储 alpha 值
            q -= alpha * y;
        }
        // 计算 Hessian 近似
        Eigen::VectorXd r = q;
        for (int i = 0; i < s_list.size(); ++i)
        {
            Eigen::VectorXd s = s_list[i];
            Eigen::VectorXd y = y_list[i];
            double rho = 1.0 / y.dot(s);
            double beta = rho * y.dot(r);
            r += s * (alpha_list[s_list.size() - 1 - i] - beta);
        }
        // 计算搜索方向 p = -r
        Eigen::VectorXd p = -r;
        // 线搜索：固定步长为1.0，可以根据需要实现更复杂的线搜索方法
        double alpha = 1.0;
        Eigen::VectorXd X_new = X + alpha * p;
        // 计算新的梯度
        Eigen::VectorXd grad_new = gradient(A, X_new, B);
        // 计算差值 s = X_new - X 和 y = grad_new - grad
        Eigen::VectorXd s = X_new - X;
        Eigen::VectorXd y = grad_new - grad;
        // 更新差分对
        if (s_list.size() == m)
        {
            s_list.erase(s_list.begin());
            y_list.erase(y_list.begin());
        }
        s_list.push_back(s);
        y_list.push_back(y);
        // 更新 X 和 grad
        X = X_new;
        grad = grad_new;
        iter++;
    }
    // 输出结果
    std::cout << "解 X: " << X.transpose() << std::endl;
}
avatar
文章
阅读
粉丝