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

基本上我需要上述。我已经浏览了谷歌，并找不到实现这一点的方法。复对称三对角矩阵的快速矩阵指数

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

让我晕眩的是Mathematica能够快速计算矩阵指数（无论它是否矩阵tridiag我不知道）。

有人能帮忙吗？

编辑：这是我想出来的，有什么意见？希望对未来的读者有用

我一直在C++的一段时间，所以下面的代码可能有点scrappy /慢，所以请赐教，如果你看到改进。

//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; 
}

来源

2012-04-12 Quantum_Oli

你需要一个指数还是多个乘法器？在后一种情况下，建议将矩阵对角线化，Mathematica可能会自动进行此操作。 – leftaroundabout 2012-04-12 13:49:30

许多不同的乘法。我有一个实际的三对角矩阵U，并且需要计算expm（i z U），其中对于许多z值，其中i = sqrt（-1）。如果你推荐对角线化矩阵，你可以推荐一个快速的C++库来做到这一点吗？ – 2012-04-12 15:01:34

这就是我的想法。我个人使用CULA来处理这个问题，这是专门针对Nvidia GPU的。 GSL也有不错的算法，应该工作得很好。然而，两者都是相当低级的 - 没有很好的C++接口。 – leftaroundabout 2012-04-12 15:20:10

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

来源

2012-04-16 11:34:26

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

回答

相关问题