c++ - Armadillo C++ expmat 堆栈

Question

我正在尝试执行 Armadillo c++ 库中的 expmat() 函数，但是当我运行它时它会被堆叠起来。首先，我在 Matlab 中尝试使用 expm 函数，并在 C++ 程序中使用不同的值，它在两种情况下都有效。所以，问题出在与我放在矩阵上的值相关的东西上，似乎其中一些太小了......这是我想要计算指数的 [14x14] 平方矩阵:

mat A;
A  << 0 << -0.0000769006 << -0.0000511322 << -0.0000915495 << 0 << 0 << -0.00254523 <<                 0.00000001 << 0 << 0 << 0 << 0 << 0 << 0 << endr
    << 0.0000769006 << 0 << 0.0000915495 << -0.0000511322 << 0.0043037 << -0.00254523 << 0 << 0 << 0.00000001 << 0 << 0 << 0 << 0 << 0 << endr
    << 0.0000511322 << -0.0000915495 << 0 << 0.0000769006 << 0.00254523 << 0.0043037 << 0 << 0 << 0 << 0.00000001 << 0 << 0 << 0 << 0 << endr
    << 0.0000915495 << 0.0000511322 << -7.69006/100000 << 0 << 0 << 0 << 0.0043037 << 0 << 0 << 0 << 0.00000001 << 0 << 0 << 0 << endr
    << 0 << 0 << 0 << 0 << 0 << 0 << 0 << 0 << 0 << 0 << 0 << 0.00000001 << 0 << 0 << endr
    << 0 << 0 << 0 << 0 << 0 << 0 << 0 << 0 << 0 << 0 << 0 << 0 << 0.00000001 << 0 << endr
    << 0 << 0 << 0 << 0 << 0 << 0 << 0 << 0 << 0 << 0 << 0 << 0 << 0 << 0.00000001 << endr
    << 0 << 0 << 0 << 0 << 0 << 0 << 0 << 0 << -0.0000769006 << -0.0000511322 << -0.0000915495 << 0 << 0 << 0 << endr
    << 0 << 0 << 0 << 0 << 0 << 0 << 0 << 0.0000769006 << 0 << 0.0000915495 << -0.0000511322 << 0 << 0 << 0 << endr
    << 0 << 0 << 0 << 0 << 0 << 0 << 0 << 0.0000511322 << -0.0000915495<< 0 << 0.0000769006 << 0 << 0 << 0 << endr
    << 0 << 0 << 0 << 0 << 0 << 0 << 0 << 0.0000915495 << 0.0000511322 << -0.0000769006 << 0 << 0 << 0 << 0 << endr
    << 0 << 0 << 0 << 0 << 0 << 0 << 0 << 0 << -0.0043037 << -0.00254523 << 0 << 0 << 0 << 0 << endr
    << 0 << 0 << 0 << 0 << 0 << 0 << 0 << 0 << 0.00254523 << -0.0043037 << 0 << 0 << 0 << 0 << endr
    << 0 << 0 << 0 << 0 << 0 << 0 << 0 << 0.00254523 << 0 << 0 << -0.0043037 << 0 << 0 << 0 << endr;

B = expmat(A);

当我在调试中按下暂停键时，我看到正在执行的行是:

arma_fortran(arma_dgemm)(transA, transB, m, n, k, (const T*)alpha, (const T*)A, ldA, (const T*)B, ldB, (const T*)beta, (T*)C, ldC);  

...在文件“blas_wrapper.hpp”中。没有显示错误，只有程序在那里，但它不会前进。如果我在矩阵上使用不同的值，我可以让它工作，但是我显示的那些值不......

关于LAPACK版本我不知道怎么查...

谢谢，

AV9

Answer 1

事实上，任何具有小于 1 的无限范数的矩阵都可能触发某种近乎无限的循环。问题不是来自您的代码，而是来自库。

以下由 g++ main.cpp -o main -larmadillo 编译的代码重现了这个问题:

#include <iostream>
#include <cstdio>
#include <armadillo>

using namespace std;
using namespace arma;

int
main( int argc, char* argv[] )
{
    cout<<"start"<<endl;

    mat A=randu(2,2);
    A=A/100;
    cout<<"A is "<<endl<<A<<endl;

    double norm_val = norm(A, "inf");  
    double log2_val = eop_aux::log2(norm_val); 
    const uword s = (std::max)(uword(0), uword(log2_val) + uword(1) + uword(1));

    cout<<"norm inf is "<<norm_val<<" log of norm inf is "<<log2_val <<" uword "<<uword(log2_val)<<" nb it is "<<s<<endl;

    mat B = expmat(A);

    cout<<"B is "<<endl<<B<<endl;
    cout<<"end"<<endl;
}

这就是为什么......

Armadillo 引用了以下文章来计算矩阵的指数:Nineteen Dubious Ways to Compute the Exponential of a Matrix, Twenty-Five Years Later由 C. Moler 和 C. Van Loan 在 SIAM Review 中发表，2003 年。第三种方法称为“缩放和平方”，由 Armadillo 图书馆实现。使用以下公式:

m 应该足够大，以便降低 A/m 的范数并缓解 A/m 的泰勒级数指数化.

算法如下:

• 选择整数s使得A/2^s的无限范数小于1
• mat AA=A/2^s
• 通过泰勒级数计算e^(AA)
• 通过对 e^(AA) s 次求平方来计算 e^(A)

在库 Armadillo 的文件 armadillo-4.550.3/include/armadillo_bits/op_expmat_meat.hpp 中，s 的计算方式如下:

    double norm_val = norm(A, "inf");  
    double log2_val = eop_aux::log2(norm_val); 
    const uword s = (std::max)(uword(0), uword(log2_val) + uword(1) + uword(1));

如果 A 的无限范数小于 1，则 log2_val 为负，将其转换为 uword(u 代表无符号)失败。它导致 s 非常大(例如 4294967291)，并且在执行平方时代码陷入几乎无限循环。

应将测试添加到库中，在文件 armadillo-4.550.3/include/armadillo_bits/op_expmat_meat.hpp 中，类似于:

uword s=42;
if(log2_val<0){
   //no need to use the formula, use Taylor series
   s=0;
 }else{
   //need to decrease A and use the formula
   s=uword(log2_val)+2;
 }

代替:

const uword s = (std::max)(uword(0), uword(log2_val) + uword(1) + uword(1));

更改文件 op_expmat_meat.hpp，通过键入 cmake 重新安装 Armadillo 。 make 和 make install(或 sudo make install) 问题就会消失。新的结果似乎是正确的。如果 A 接近 0:

我将通过电子邮件将此错误报告给 Ar​​madillo 的开发人员。编辑:不需要这样做:正如@mtall 所注意到的，此问题已在版本 4.550.4 中删除。我的猜测是开发人员已经看到了您的帖子。几个小时前修改了补丁版本:开发人员的团队 react 迅速!