复对称三对角矩阵的快速矩阵指数答案

【问题标题】：Fast matrix exponential of a complex symmetric tridiagonal matrix复对称三对角矩阵的快速矩阵指数
【发布时间】：2012-04-12 13:38:24
【问题描述】：

基本上我需要上面的。我已经搜索了谷歌，但找不到实现这一点的方法。

我在这里http://www.guwi17.de/ublas/examples/ 找到了这个功能，但它太慢了。我什至按照 MATLAB 的例程编写了自己的 Pade Approximation，但它只比链接中的快一点。

让我震惊的是 Mathematica 计算矩阵指数的速度有多快（我不知道它是否关心矩阵是否是 tridiag）。

有人可以帮忙吗？

编辑：这是我想出的，任何 cmets？希望对未来的读者有用

我已经使用 c++ 一段时间了，所以下面的代码可能有点杂乱/缓慢，所以如果你看到改进，请赐教。

//Program will compute the matrix exponential expm( i gA ). where i = sqrt(-1) and gA is a matrix

#include <iostream>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>


int main() {

    int n = 100;

    unsigned i = 0;
    unsigned j = 0;

    gsl_complex img = gsl_complex_rect (0.,1.);

    gsl_matrix * gA = gsl_matrix_alloc (n, n);


//Set gA:       
    for ( i = 0; i<n; i++) {
        for ( j = 0; j<=i; j++) {
            if( i == j ) {
                if( i == 0 || i == n-1 ) {
                    gsl_matrix_set (gA, i, j, 1.);
                } else {
                    gsl_matrix_set (gA, i, j, 2.);
                }
            } else if( j == i-1 ) {
                gsl_matrix_set (gA, i, j, 1.);
                gsl_matrix_set (gA, j, i, 1.);
            } else {
                gsl_matrix_set (gA, i, j, 0.);
                gsl_matrix_set (gA, j, i, 0.);
        }
    }
}


//get SVD of gA 
    gsl_matrix * gA_t = gsl_matrix_alloc (n, n);
    gsl_matrix_transpose_memcpy (gA_t, gA);                 

    gsl_matrix * U = gsl_matrix_alloc (n, n);
    gsl_matrix * V= gsl_matrix_alloc (n, n);
    gsl_vector * S = gsl_vector_alloc (n);


    gsl_vector * work = gsl_vector_alloc (n);
    gsl_linalg_SV_decomp (gA_t, V, S, work);
    gsl_vector_free(work);

    gsl_matrix_memcpy (U, gA_t);


//Take exponential of S and form matrix
    gsl_matrix_complex * Sp = gsl_matrix_complex_alloc (n, n);
    gsl_matrix_complex_set_zero (Sp);
    for (i = 0; i < n; i++) {
        gsl_matrix_complex_set (Sp, i, i, gsl_complex_exp ( gsl_complex_mul_real ( img, gsl_vector_get(S, i) ) ) ); // Vector 'S' to matrix 'Sp'
    }

    gsl_matrix_complex * Uc = gsl_matrix_complex_alloc (n, n);


//convert U to a complex matrix for next step   
    for (i = 0; i < n; i++) {
        for ( j = 0; j<n; j++) {
            gsl_matrix_complex_set (Uc, i, j, gsl_complex_rect ( gsl_matrix_get(U, i, j), 0. ) );       
        }
    }


//recombine U S Utranspose  
    gsl_matrix_complex * tA = gsl_matrix_complex_alloc (n, n);
    gsl_matrix_complex * eA = gsl_matrix_complex_alloc (n, n);
    gsl_blas_zgemm (CblasNoTrans, CblasNoTrans, GSL_COMPLEX_ONE, Uc, Sp, GSL_COMPLEX_ZERO, tA);
    gsl_blas_zgemm (CblasNoTrans, CblasTrans, GSL_COMPLEX_ONE, tA, Uc, GSL_COMPLEX_ZERO, eA);


//Print result  
    std::cout << "eA:" << std::endl;
    for (i = 0; i < n; i++)  
        for (j = 0; j < n; j++)
            printf ("%g + %g i\n", GSL_REAL(gsl_matrix_complex_get (eA, i, j)), GSL_IMAG(gsl_matrix_complex_get (eA, i, j)));
    std::cout << "\n" << std::endl;


//Free up
    gsl_matrix_free(gA_t);
    gsl_matrix_free(U);
    gsl_matrix_free(gA);
    gsl_matrix_free(V);
    gsl_vector_free(S);
    gsl_matrix_complex_free(Uc);
    gsl_matrix_complex_free(tA);    
    gsl_matrix_complex_free(eA);    

    return 0;
}

【问题讨论】：

您只需要一个指数，还是需要多个具有不同乘数的指数？在后一种情况下，建议对矩阵进行对角化，Mathematica 可能会自动执行此操作。
许多具有不同的乘数。我有一个真正的三对角矩阵 U 并且需要计算 expm(i z U) 其中 i = sqrt(-1) 的许多 z 值。如果你建议对矩阵进行对角化，你能推荐一个快速的 c++ 库吗？
我就是这么想的。我个人为此使用 CULA，它专门用于 Nvidia GPU。 GSL 也有不错的算法，应该可以很好地工作。然而，两者都是相当低级的——没有好的 C++ 接口。
我已经对角化了！往上看。它比使用我的链接中给出的 Pade 方法的函数快大约 100 倍。对于我的手段，这应该足够快。我也只需要进行一次分解。感谢您的帮助。

标签： c++ c matrix blas

【解决方案1】：

上面我发布到我的问题中的代码效果很好。我发现的另一个改进是通过特征值缩放 V 副本中的列，而不是执行完整的矩阵乘法。由于 zgemm 是此代码中最慢的部分，因此删除其中一个 zgemm 函数将其速度提高了近两倍。我很快就会在这里发布完整的代码清单。

【讨论】：