Armadillo C ++ expmat堆栈 [英] Armadillo C++ expmat stacks

问题描述

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

I'm trying to execute the expmat() function from the Armadillo c++ library but it gets stacked when i run it. First i tried in Matlab with the expm function, and with different values in the C++ program and it works in both cases. So, the problem is with something related with the values i put on the matrix, seem that some of them are too small... This is the squared [14x14] matrix i want to calculate the exponential:

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);

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

When i press pause on the debug i see that line wich is being executed is:

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"中.不显示任何错误，仅显示程序在其中，但不会前进.如果我在矩阵上使用不同的值，则可以使其正常工作，但是显示的值不正确...

...inside the file "blas_wrapper.hpp". No errors are displayed, only the program is there but it doesn´t advance. If i use different values on the matrix i can make it working, but with the ones i displayed not...

关于LAPACK版本，我不知道如何检查...

About the LAPACK version i don´t know how to check it...

谢谢

AV9

推荐答案

实际上，任何无限范数小于1的矩阵都可能触发某种几乎无限的循环.问题并非来自于您的代码，但来自库.

In fact, any matrix with infinite norm lower than 1 may trigger some sort of nearly infinite loop. The problem does not come from your code, but from the library.

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

The following code, compiled by g++ main.cpp -o main -larmadillo reproduces the problem :

#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引用了以下文章来计算矩阵的指数:十九种可疑的方法C. Moler和C. Van Loan在2003年的SIAM评论中计算矩阵的指数，二十五年后..Armadillo库实现了第三种方法，称为缩放和平方".使用以下公式:

Armadillo quotes the following article to compute the matrix's exponential : Nineteen Dubious Ways to Compute the Exponential of a Matrix, Twenty-Five Years Later by C. Moler and C. Van Loan in SIAM Review, 2003. The third method, called "Scaling and squaring" is implemented by the library Armadillo. The following formula is used :

m应该足够大，以免破坏A/m的范数并简化泰勒级数对A/m的指数化.

m should be large enough so as to decrase the norm of A/m and ease the exponentialisation of A/m by Taylor series.

算法如下:

• 选择整数s，以使A/2^s的无穷范数小于1.
• mat AA=A/2^s
• 通过泰勒级数计算e^(AA)
• 通过计算e^(AA) s次来计算e^(A)

• Choose integer s so that the infinite norm of A/2^s is lower than one
• mat AA=A/2^s
• compute e^(AA) by Taylor series
• compute e^(A) by squarring e^(AA) s times

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

In the library Armadillo, in file armadillo-4.550.3/include/armadillo_bits/op_expmat_meat.hpp here is how s is computed :

    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)，并且在执行平方运算时，代码陷入了几乎无限的循环中.

If the infinite norm of A is lower than 1, log2_val is negative and casting it to uword (u stands for unsigned) fails. It results in s being very large (4294967291 for instance), and the code is stuck in a nearly infinite loop when squaring is performed.

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

A test should be added to the library, in file armadillo-4.550.3/include/armadillo_bits/op_expmat_meat.hpp, something like :

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 . make和make install(或sudo make install)重新安装犰狳，该问题将消失.新的结果似乎是正确的.如果A接近0:

Change file op_expmat_meat.hpp, reinstall armadillo by typing cmake . make and make install (or sudo make install) and the issue will be gone. The new result seems to be correct. If A is near 0 :

我将通过电子邮件向Armadillo的开发人员报告此错误. 无需这样做:@mtall注意到，此问题在4.550.4版中已删除.我的猜测是开发人员已经看过您的帖子.补丁版本已在数小时前进行了修改:开发人员团队处于反应状态！

I will report this bug to the developpers of Armadillo by email. EDIT : no need to do so : as noticed by @mtall, this issue is removed in version 4.550.4. My guess is that the developpers have seen your post. The patched version has been modified a few hours ago : the developpers'team is reactive !

这篇关于Armadillo C ++ expmat堆栈的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！