2
$\begingroup$

I would like to evaluate a polynomial of matrix $g(A)= \sum_{m=0}^{k}e^{-i m \phi} A^{m}$ where $\phi$ is some angle provided by the user and $k$ is some positve integer. So I did a very naive thing to evaluate it

gkm[k_,\[Phi]_,x_]:= Module[{mpow}, mpow= IdentityMatrix[Length[x]]; result= mpow; Table[ mpow= mpow . x; result+= Exp[-I m \[Phi]]mpow;, {m,1,k}]; Return[result] ]

Basically, I am saving the result in each iteration so that I do not have to multiply $m$ times to get $U^{m}$ . However, this is not really optimal. So for eg. if $x$ were not matrices then the polynomial is efficiently evaluated using Horner's method. However, as mentioned in Wikipedia, it is not optimal if $x$ is a matrix. From there I went to the method of preprocessing, however the method completely eludes me. I was wondering if there is an efficient way to evaluate polynomial $g(A)$ faster than what I already have. I require it for matrices of dimensions $2^{16}$ and $k \sim 12$ ( which I would run on the workstation that I have access to).

Edit 1: I understand that the series is a geometric series, so it can be exactly summed as mentioned in the comment below. However, I would like to avoid finding an inverse at any cost. Another way that I already tried is using the MatrixFunction command of Mathematica and weirdly, it gives an incorrect result even for matrices of dimension $2^{5}$ (I have no idea why). I am assuming it is doing some kind of approximation scheme, but it is a bad one, it seems (atleast to me), though it is surely faster than what I already have.

As a benchmark, the code that I wrote already does the job for a $2^{10}$ dimensional matrix with $k=4$ in $\sim 7.5s$ (Macbook Air M1). Matrix $x$ is fully dense, complex with floating point entries.

Edit2: I tried to use LinearSolve for the case when I needed to calculate the inverse (closed result for the sum). It is not bad, with slightly more memory requirements than the gkm function that I wrote. It seems that I need to reduce the number of matrix-matrix multiplications as possible. So maybe Horner's method would be a way to go at least. Here what I tried and might give a little bit better than what I already have. If a polynomial $p(x)$ is given then Mathematica can already give the Horner's form using HornerForm. However, now I am a bit stuck on how to use the resulting form for matrices. As an example consider $k=4, \phi= \pi/2$ then I can do the following

gk[k_,\[Phi]_,x_]:=Sum[Exp[-I m \[Phi]] x^m,{m,0,k}] HornerForm[gk[k, \[Pi]/2, x]] 1+x (-I+x (-1+x (I+x)))

But now I want $x$ to be a matrix and $1$ to be the identity matrix, and the simple product should be changed to dot product for matrices. Is there a way to convert it into a form such that it can be used for matrices? Yes there is and I found the way.

 (HornerForm[gk[k,\[Pi]/2,y]]/.{Times-> Dot,-I -> -I IdentityMatrix[2^\[CapitalNu]],I->I IdentityMatrix[2^\[CapitalNu]],1-> IdentityMatrix[2^\[CapitalNu]],-1-> -IdentityMatrix[2^\[CapitalNu]],y->U})

The above substitution takes care of it, and I have checked it. The time taken is similar to the first code that I wrote. so it seems it is still the same.

Improve this question
$\endgroup$
20
  • 1
    $\begingroup$ First of all, your exponential prefactor can be incorporated in the definition of $A$. Thus, you question is reduced to a geometric series: $S_k=\sum_m^{k}A^m=(I-A)^{-1}(I-A^{k+1})$. This holds provided the matrix is invertible. $\endgroup$ Commented Mar 2 at 9:36
  • $\begingroup$ Agreed, but it is not really "efficient." to invert a matrix of large dimensions. Please feel free to correct me. I am dealing with matrices as large as $2^{16}$.) $\endgroup$ Commented Mar 2 at 10:37
  • 1
    $\begingroup$ How are you dealing with matrices of the size $2^{16}$ on a Macbook Air M1 if it ha maximal configurable memory of 16GB and a real matrix (8 byte per element) of the size $2^{16}$ requires 32GB? You need to keep the original matrix and the result, which means you need at leas 64 GB. Please, specify the problem realistically. $\endgroup$ Commented Mar 2 at 13:45
  • 1
    $\begingroup$ The answer may depend on what type of entries in the matrix x you have. Are those exact numbers or floating point ones? $\endgroup$ Commented Mar 3 at 1:47
  • 1
    $\begingroup$ It is likely that the reason for MatrixFunction's poor accuracy here is the infamous sign problem, for which AFAIK the only general remedy is more accuracy. MatrixFunction is probably using actual floats, whereas Mathematica will in other contexts generally do things exactly if it can and can be told to use arithmetic to some specific precision otherwise. $\endgroup$ Commented Mar 3 at 2:21

3 Answers 3

Reset to default
3
$\begingroup$

As indicated by many comments, optimal answer crucially depends on the size of input data --- tradeoff between computational time and memory efficiency.

For moderate size of matrices, I think the build in function MatrixPower is reasonably fast:

ClearAll[A, k, n]; n = 2^10; k = 12; phi = Pi/7; A = RandomReal[{-10, 10}, {n, n}]; B = Exp[-I*phi]*A; AbsoluteTiming[Sum[MatrixPower[B, m], {m, 0, k}];] (* {1.28058, Null} *) $Version (* "14.1.0 for Mac OS X ARM (64-bit) (July 16, 2024)" *)
Improve this answer
$\endgroup$
1
$\begingroup$

Note that MatrixPower[A,m] will likely be quite densely populated if n us large enough, even if A is quite sparse. (The set matrices for which this does not happen is relatively small.

Hence MatrixPower[A,m] will likely have complexity $O(n^3)$ for an $n \times n$ matrix A and if m is large enough.

Often you need not $g(A)$ itself, but rather its action $v \mapsto g(A) \, v$ on vectors. This routine does this for you: DotA is a function that implements the operation of A on vectors, coefficients is the coefficient array of g and v is a vector.

MatrixPolynomial[DotA_Function, coefficients_?VectorQ, v_?VectorQ] := Module[{degree, Akv, result}, degree = Length[coefficients] - 1; Akv = v; result = coefficients[[1]] v; Do[ Akv = DotA[Akv]; result += coefficients[[k + 1]] Akv; , {k, 1, degree}]; result ]

Here a usage example featuring a sparse matrix A.

A = KirchhoffMatrix[RandomGraph[{2^12, 2^15}]]; DotA = v |-> A . v; g = t |-> 1 + t + t^2 + t^3 + t^4 + t^5 + t^6 + t^7 + t^8 + t^9 + t^10; rules = CoefficientRules[g[t], t]; coefficients = Partition[Range[0, Max[Keys[rules]]], 1] /. rules; degree = Length[coefficients] - 1;

Using dense matrix powers:

gATime = First@AbsoluteTiming[ Block[{A = Normal[A], B}, B = IdentityMatrix[Length[A], WorkingPrecision -> MachinePrecision]; gA = coefficients[[1]] B + Sum[B = B . A; coefficients[[k + 1]] B, {k, 1, Length[coefficients] - 1}]; ] ]

4.6044

Using matrix-vector products only:

v = RandomReal[{-1, 1}, Length[A]]; gAvTime = First@RepeatedTiming[ gAv = gA . v; ] gAvMatriFreeTime = First@RepeatedTiming[ gAvMatriFree = MatrixPolynomial[DotA, coefficients, v]; ] relativeError = Norm[gAv - gAvMatriFree]/Norm[gAv]

0.0166871

0.000680698

8.64671*10^-16

Plot[{gATime + m gAvTime, m gAvMatriFreeTime}, {m, 1, 100}, PlotLabel -> "Time for m matrix-vector multiplications with g[A]", AxesLabel -> {"m", "t"}, PlotLegends -> {"dense matrix", "matrix free"} ]

Time for m matrix-vector multiplications with g[A]

Improve this answer
$\endgroup$
1
$\begingroup$

For any complex x we have $$\sum _{m=0}^{\infty } x^m e^{i m \phi }-\sum _{m=k}^{\infty } x^m \ e^{i m \phi } = \frac{1-x^{k} e^{ i \ k \ \phi }}{1-x \ e^{i \ \phi }}$$

Since all matrices commute, define

 f[x_, \[Phi]_, m_] := ((# - MatrixPower[E^(I \[Phi]) x, m]) . MatrixPower[# - E^(I \[Phi]) x, -1] &) [IdentityMatrix[Length[x]]] A = RandomComplex[{-1 - I, 1 + I}, {1024, 1024}]; Timing[f[A, 0.98 \[Pi], 12];] {1.29688, Null} MaxMemoryUsed[]-MemoryInUse[] 222 MB

on a notebook Core I5 10th gen

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.