Trigonometry

Overview

In the trigonometric solvers module, we provide support for computing the sine and cosine of matrices.

Method

The basis of these methods is to expand the cosine function on a set of Chebyshev polynomials, and use that to compute the matrix function. Using the base cosine function, it is of course possible to also compute the sine function, using the identity:

One difficulty of this approach is that the spectrum of the input matrix not necessarily contained in the range [-1, 1]. This can be addressed using the double angle formula:

The cosine function is an even function, which allows for further savings when expanding the matrix using Chebyshev polynomials.

Example Implementation

def DivideChebCompute(InputMatrix):
  # Determine the Scaling Factor
  matrix_dimension = InputMatrix.shape[0]
  min_val, max_val = Helper.EigenCircle(InputMatrix)
  spectral_radius = max_val
  sigma_val = 1
  sigma_counter = 1
  while spectral_radius/sigma_val > 1.0:
    sigma_val = sigma_val*2
    sigma_counter = sigma_counter + 1
  ScaledMatrix = InputMatrix/sigma_val

  # Compute the chebyshev coefficients
  x = linspace(-1,1,2000)
  y = [cos(i) for i in x]
  coefficient_list = chebfit(x,y,17)

  identity_matrix = identity(matrix_dimension)
  # Basic
  T0 = scipy.sparse.identity(matrix_dimension)
  T1 = ScaledMatrix.copy()
  T2 = 2.0*ScaledMatrix.dot(T1) - T0

  # First Half
  T4 = 2.0*T2.dot(T2) - T0
  T6 = 2.0*T4.dot(T2) - T2
  T8 = 2.0*T6.dot(T2) - T4

  # Compute the second half
  OutputMatrix = 0.5*coefficient_list[16]*T8
  OutputMatrix = OutputMatrix + 0.5*coefficient_list[14]*T6
  OutputMatrix = OutputMatrix + 0.5*coefficient_list[12]*T4
  OutputMatrix = OutputMatrix + 0.5*coefficient_list[10]*T2
  OutputMatrix = T8.dot(OutputMatrix)

  # Add in the first half
  OutputMatrix = OutputMatrix + (coefficient_list[8])*T8
  OutputMatrix = OutputMatrix + (coefficient_list[6]+coefficient_list[10]/2.0)*T6
  OutputMatrix = OutputMatrix + (coefficient_list[4]+coefficient_list[12]/2.0)*T4
  OutputMatrix = OutputMatrix + (coefficient_list[2]+coefficient_list[14]/2.0)*T2
  OutputMatrix = OutputMatrix + (coefficient_list[0]+coefficient_list[16]/2.0)*T0

  for i in range(0,sigma_counter-1):
    OutputMatrix = 2.0*OutputMatrix.dot(OutputMatrix) - identity_matrix
  return OutputMatrix

How To Cite

It's appropriate to cite the following two articles for the basic double angle method:

@article{serbin1980algorithm,
title={An algorithm for computing the matrix cosine},
author={Serbin, Steven M and Blalock, Sybil A},
journal={SIAM Journal on Scientific and Statistical Computing},
volume={1},
number={2},
pages={198--204},
year={1980},
publisher={SIAM}
}

@article{higham2003computing,
title={Computing the matrix cosine},
author={Higham, Nicholas J and Smith, Matthew I},
journal={Numerical Algorithms},
volume={34},
number={1},
pages={13--26},
year={2003},
publisher={Springer}
}

The following article was the inspiration for the particular factorization utilized:

@article{yau1993reducing,
title={Reducing the symmetric matrix eigenvalue problem to matrix multiplications},
author={Yau, Shing-Tung and Lu, Ya Yan},
journal={SIAM Journal on Scientific Computing},
volume={14},
number={1},
pages={121--136},
year={1993},
publisher={SIAM}
}