# Computing the Matrix Logarithm with the Romberg Integration Method
## 1 Introduction
The inverse scaling and squaring procedure, initially proposed in [27], is an extension to the matrix domain of the technique used by [PERSON] to compute his table of logarithms, collected in [45]. This method takes advantage of the matrix identity \(\log\left(A\right)=2^{s}\log\left(A^{2^{-s}}\right)\) and evaluates \(\log\left(A\right)\) by combining argument reduction and approximation. By taking a certain number of square roots of \(A\), the problem is reduced to the computation of the logarithm of a matrix with eigenvalues close to \(1\). Indeed, the approximation of the matrix logarithm is performed in three different stages, as explained in [27]:  
1. Find an integer \(s\), so that matrix \(A^{2^{-s}}\) is close to the identity matrix \(I\). For that purpose, an algorithm that computes matrix square roots must be employed;
2. Approximate \(\log\left(A^{2^{-s}}\right)\) by \(r_{m}(A^{2^{-s}}-I)\), where \(r_{m}\) is the diagonal Pade approximant of degree \(m\) to the function \(\log\left(1+x\right)\);
3. Compute the approximation \(\log\left(A\right)\approx 2^{s}r_{m}(A^{2^{-s}}-I)\). Taking into account the previous three-stage procedure and the following integral expression for the logarithm [40]  
\[\log(A)=\int_{0}^{1}(A-I)((A-I)x+I)^{-1}dx, \tag{2}\]  
our proposal in this work is to compute the matrix logarithm in three phases somewhat similar to those described above. Notwithstanding, \(\log\left(A^{2^{-s}}\right)\) is approximated in the second stage by means of the expression (2) and using the well-known Romberg method [46]. In addition, the matrix square roots in the first phase are worked out thanks to the scaled Denman-Beavers iteration explained in [1]. The importance of the method used for matrix square root computation must be emphasized.Working with one method or another has a direct impact, not only on the accuracy of the result, but also on the associated execution time. [PERSON]'s method employs the trapezoidal rule to approximate numerically the definite integral \(I(f)=\int_{a}^{b}f(x)dx\). It reduces the integration step \(h\) by half at each iteration and applies the Richardson extrapolation formula to the previous results. This quadrature method generates the so-called Romberg tableau, in the form of a lower triangular matrix \(R\), whose elements are numerical estimates of the definite integral to be calculated. If we take  
\[h_{i}=\frac{b-a}{2^{i-1}},\quad i\geq 1, \tag{3}\]as the integration step to be applied in the successive iterations, we define the calculation of each term \(\{R_{i,j}\}\) as  
\[\left. \begin{array}{rcl}R_{1,1}&=&\frac{b-a}{2}(f(a)+f(b))\\ R_{i,1}&=&\frac{R_{i-1,1}}{2}+h_{i}\sum_{k=1}^{2i-2}f(a+(2k-1)h_{i}),i\geq 2\\ R_{i,j}&=&\frac{4^{j-1}R_{i,j-1}-R_{i-1,j-1}}{4^{j-1}-1},2\leq i,2\leq j\leq i \end{array} \right\}. \tag{4}\]  
Below is a diagram depicting how the triangular matrix elements are filled after \(m\) iterations, where the arrows indicate the order in which each term is generated from the previous ones:  
\[\left. \begin{array}{ccccccccc}R_{1,1}&&&&&&&&\\ \downarrow&\searrow&&&&\\ R_{2,1}&\rightarrow&R_{2,2}&&&&\\ \downarrow&\searrow&&\searrow&&\\ R_{3,1}&\rightarrow&R_{3,2}&\rightarrow&R_{3,3}&&&&\\ \downarrow&\searrow&&\searrow&&\searrow&&\\ R_{4,1}&\rightarrow&R_{4,2}&\rightarrow&R_{4,3}&\rightarrow&R_{4,A}&&\\ \vdots&&\vdots&&\vdots&&\vdots&\ddots&\\ R_{m,1}&\rightarrow&R_{m,2}&\rightarrow&R_{m,3}&\rightarrow&R_{m,A}& \ldots&R_{m,m}.\end{array} \right. \tag{5}\]  
The last diagonal term \(R_{m,m}\) of the matrix is always the most accurate estimate of the integral. In addition, according to ([46] p. 342), it is known that if function \(f\) has \(2m+2\) continuous derivatives, then the asymptotic error is  
\[|I(f)-R_{m,m}|=O\Big{(}4^{-m}h^{2m+2}\Big{)}. \tag{6}\]  
Moreover, if \(m\) is sufficiently large, then  
\[|I(f)-R_{m,m}|<|R_{m,m}-R_{m-1,m-1}|. \tag{7}\]  
Therefore, the previous process must be continued until the difference between two successive diagonal elements becomes sufficiently small.