# Computing the Matrix Logarithm with the Romberg Integration Method
## 2 Romberg Method for the Matrix Logarithm
### Theoretical Analysis of the Error
19:\(L=2^{2}L\)\(\triangleright\) Squaring the logarithm matrix
20:\(L=TLT^{-1}\)\(\triangleright\) Postprocessing of matrix \(L\)
```  
**Algorithm 1** Given a matrix \(A\in\mathbb{C}^{n\times n}\), a maximum number of square roots \(max\_sqrts\) to be computed, a maximum value \(m\in\{1,2,3,\cdots,15\}\), and a tolerance \(tol\), this algorithm computes \(L=\log(A)\) by the Romberg approximation and the inverse scaling and squaring technique. As an aside, it should be noted that it is feasible to estimate the \(1\)-norm of \((A-I)^{2m+1}\) in the Steps 3 and 6 more quickly if we understand that  
\[E_{2}=I(f)-R_{2,2}=\sum\limits_{k\geq 2}a_{k}^{(2)},\]0  
taking into account that \(A^{2}-2A+I\) can be evaluated without matrix products, as long as at least one square root of \(A\) has been previously computed (and saving matrix \(A\) before computing the square root). A similar improvement can be carried out in the Step 12. Algorithm 2 works out the matrix logarithm using the Romberg method. In Steps 1 to 5, it initializes the suitable variables and provides the first diagonal term in \(R_{1}\), following the expression in (9). At each iteration \(i\), within the while-loop that involves Steps 6 to 21, the integration step \(h\) is halved, and the \(i\)-th row of the Romberg tableau is completed in matrices \(R_{1}\) to \(R_{i}\). It should be noted that the algorithm does not keep in memory all the elements that would compose the corresponding whole block of the lower triangular matrix, in a attempt to optimize computer memory and execution time. To achieve this purpose, each new element of each row \(i\) is saved in the auxiliary matrices \(T_{1}\) and \(T_{2}\) alternately, according to Formulas (9) and (11), and matrices \(R_{1}\) to \(R_{i-1}\) are reused to store the terms of the row currently in process, once the components of the previous one are no longer needed. In Step 18, the convergence criterion is evaluated by comparing the diagonal term \(R_{i}\) with the equivalent element of the preceding row, as previously preserved in matrix \(L\). The storage and computational costs of the two proposed algorithms, together with the codes on which they depend, are given in \[E_{2}=I(f)-R_{2,2}=\sum\limits_{k\geq 2}a_{k}^{(2)},\]1 The cost of storage is expressed as the maximum number of matrices required to be stored together in memory. The computational cost is measured in terms of the number of matrix products. Other lineal algebra operationssuch as systems of linear equations and matrix inverses also appear in the expression for the cost. It is assumed that either of these two operations has a cost of \(4/3\) matrix products [50]. ```
1:\(h=1\)
2:\(R_{1}=0.5(f(0,A)+f(1,A))\)\(\triangleright\) Where \(f(x,A)=(A-I)((A-I)x+I)^{-1}\)
3:\(L=R_{1}\)
4:\(i=2\)
5:\(finish=0\)
6:while\(finish=0\)and\(i\leq m\)do
7:\(h=h/2\)
8:\(S=0\)
9:for\(k=1:2^{i-2}\)do
10:\(S=S+f((2k-1)h,A)\)
11:endfor
12:\(T_{1}=0.5R_{1}+hS\)
13:for\(j=2:i\)do
14:\(T\mod(j-1,2)+1=(4^{j-1}T\mod(j-2,2)+1-R_{j-1})/(4^{j-1}-1)\)
15:\(R_{j-1}=T\mod(j-2,2)+1\)
16:endfor
17:\(R_{i}=T\mod(i-1,2)+1\)
18:\(finish=(\|R_{i}-L\|\leq tol)\)
19:\(L=R_{i}\)
20:\(i=i+1\)
21:endwhile
```  
**Algorithm 2** Given a matrix \(A\in\mathbb{C}^{n\times n}\), this algorithm computes \(L=\log(A)\) by the Romberg method, with a maximum number of iterations \(m\) and a tolerance \(tol\). The cost of the algorithms is provided by themselves or together with the various codes they invoke. For example, Algorithm 1 needs only three matrices to be stored by itself but, globally, when it is being executed and the different functions in which it is structured are called, it saves \(10+m\) matrices simultaneously in memory, at most. Regarding the computational cost, it should be noted that most of the algorithms are iterative. As an example, the Denman-Beavers code cost is expressed in terms of number of iterations \(k\) required to reach convergence. In each iteration, the inverse of two matrices must be calculated, which is equivalent to a cost of \(8/3\) matrix products.