在犰狳中将矩阵除以向量

Posted 2023-02-16

【中文标题】在犰狳中将矩阵除以向量

【英文标题】：Divide Matrix by vector in armadillo

【发布时间】：2017-06-01 13:24:34

【问题描述】：

我有一个稀疏矩阵和一个向量，它们都是用犰狳构造的。现在我想将矩阵除以向量，即将矩阵的第一列除以第一个元素，依此类推。现在我可以循环执行了：

int r_num = 5000;
arma::colvec r_vec = arma::linspace(0, 1, r_num);
double dh = abs(r_vec(1)-r_vec(0));
arma::sp_cx_mat r1_matrix = arma::sp_cx_mat(r_num, r_num);
r1_matrix.diag(0).fill(0);
r1_matrix.diag(1).fill(8);
r1_matrix.diag(-1).fill(-8);
r1_matrix.diag(2).fill(-1);
r1_matrix.diag(-2).fill(1);

arma::cx_colvec divider_vec = r_vec*std::complex<double>1.0, 0.0;
divider_vec(0) = 1;

for (size_t i = 0; i < r_num; i++)

    if (i % 100 == 0)
        std::cout << i << " of " << r_num << '\n';
    r1_matrix.col(i) = r1_matrix.col(i) / divider_vec(i);

但这非常慢（尤其是在我之后会增加r_num 之后）。当使用通常的方法时

r1_matrix = r1_matrix/divider_vec;

我收到尺寸不匹配的错误 ((5000x5000) 和 (5000x1))。还有其他更快的方法吗？

【问题讨论】：

我可能遗漏了一些东西，但是，您不能只计算向量的“逆”，然后将其与矩阵相乘吗？

你的意思是，计算r1_matrix*(1/divider_vec)?

@Zouch：这会将矩阵简化为我不想要的向量。

如果你让 divider_vec 成为对角矩阵？

我不确定稀疏矩阵和犰狳的具体内容，但是使用你的逆向量作为对角矩阵的对角线应该可以完成这项工作

【参考方案1】：

除了可能由于稀疏矩阵或犰狳（我根本不知道）而导致的细节之外，我认为解决您的问题的最佳方法是计算向量的逆，将其存储为对角线的对角线矩阵（矩阵的大小），然后将矩阵乘以计算得到的矩阵。

使用 Eigen 的一些示例：

#include <Eigen/Core>
#include <iostream>

using namespace Eigen;

int main()

    Matrix<float, 3, 3> r1_matrix;
    r1_matrix <<
        2, 3, 5,
        4, 6, 10,
        6, 9, 15;

    Matrix<float, 3, 1> divider_vec;
    divider_vec << 2, 3, 5;
    divider_vec = divider_vec.cwiseInverse();

    Matrix<float, 3, 3> divider_mat = divider_vec.asDiagonal();

    r1_matrix = r1_matrix * divider_mat;

    std::cout << std::endl << r1_matrix << std::endl;

    Matrix<float, 5, 3> r2_matrix;
    r2_matrix <<
        2, 3, 5,
        4, 6, 10,
        6, 9, 15,
        8, 12, 20,
        10, 15, 25;

    Matrix<float, 3, 5> divider_mat2 = Matrix<float, 3, 5>::Zero();
    divider_mat2.block<3, 3>(0, 0) = divider_mat;

    r2_matrix = (r2_matrix * divider_mat2).block<5, 3>(0, 0);
    std::cout << std::endl << r2_matrix << std::endl;

    SparseMatrix<float> sp_matrix(5, 5);
    std::vector<Triplet<float>> data =
    
         0, 0, 2.f ,  2, 1, 6.f ,  1, 2, 10.f ,  3, 3, 14.f ,  4, 3, 21.f ,  1, 4, 11.f 
    ;
    sp_matrix.setFromTriplets(data.begin(), data.end());

    Matrix<float, 5, 1> divider_vec2;
    divider_vec2 << 2, 3, 5, 7, 11;
    divider_vec2 = divider_vec2.cwiseInverse();
    SparseMatrix<float> divider_sp_mat(5, 5);
    data.clear();

    for (int i = 0; i < 5; ++i)
    
        data.push_back(Triplet<float> i, i, divider_vec2[i]);
    

    divider_sp_mat.setFromTriplets(data.begin(), data.end());

    sp_matrix = sp_matrix * divider_sp_mat;
    std::cout << std::endl << sp_matrix << std::endl;
    return 0;

我不确定当矩阵的行数多于列数时，除了制作方阵之外是否还有其他方法。