Podnoszenie macierzy do n-tej potęgi

Jak szybko i wydajnie podnosić macierz do potęgi n-tej? Znane są przynajmniej dwie metody.

Przygotowania wstępne

Do ładnego przedstawiania macierzy przyda się pakiet mat2tex.

library(devtools)
install_github("mat2tex", "markheckmann") #Tak nie poszło
install_github("markheckmann/mat2tex") #Tak sie udało

Oczywiście na początku trzeba go załadować. Należy także dodać \usepackage{amsmath} do nagłówka YAML dokumentu.

library(pracma)
library(microbenchmark)
library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.3.5     v purrr   0.3.4
## v tibble  3.1.6     v dplyr   1.0.8
## v tidyr   1.2.0     v stringr 1.4.0
## v readr   2.1.2     v forcats 0.5.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x purrr::cross()  masks pracma::cross()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(mat2tex)

mat2tex_options(digits=3, mtype="bmatrix", matenvir=1)

Metoda dziel i rządź

Metoda dziel i rządź opiera się na rekurencyjnym algorytmie typu dziel i rządź.

powm_d = function(M, n){
  if(n==1){
    M
  } else {
    P = powm_d(M, floor(n/2))
    if(mod(n,2) == 0){
      P %*% P
    } else {
      P %*% P %*% M
    }
  }
}

Metoda wartości własnych

Metoda ta polega na wykorzystaniu wektorów i wartości własnych macierzy.

Niech \(\textbf A\) będzie kwadratową macierzą o wymiarach \(n \times n\) z \(n\) liniowo niezależnymi wektorami własnymi \(q_i\) macierzy \(\textbf A\) (gdzie \(i=1,\dots,n\)). Wtedy \(\textbf A\) można rozłożyć na czynniki jako

\[ \textbf A = \textbf {QDQ}^{-1} \]

gdzie \(\textbf Q\) jest macierzą kwadratową \(n \times n\), której \(i\)-ta kolumna jest wektorem własnym \(q_i\) macierzy \(\textbf A\), a \(\textbf D\) jest macierzą diagonalną której elementy diagonalne są odpowiadającymi wartościami własnymi \(D_{ii} = \lambda_i\).

W R wartości własne oraz wektory własne można uzyskać stosując funkcję eigen(). Poniżej mały przykład jej zastosowania.

A = matrix(c(1,-3,-2,4), nrow = 2, ncol = 2, byrow=TRUE)
xx(xm(A))

\[ \begin{bmatrix} 1.000 & -3.000 \\ -2.000 & 4.000 \\ \end{bmatrix} \]

ev = eigen(A)

Wartości własne są w zmiennej values

xx(ev$values)

\[ \begin{bmatrix} 5.372 \\ -0.372 \\ \end{bmatrix} \]

Wektory własne zaś są w zmiennej vectors.

xx(xm(ev$vectors))

\[ \begin{bmatrix} 0.566 & -0.909 \\ -0.825 & -0.416 \\ \end{bmatrix} \]

Sprawdźmy czy wszystko jest poprawne

xx("A &= QDQ^{-1}", "\\\\", xm(A), " &= ", xm(zapsmall((ev$vectors %*% diag(ev$values)) %*% inv(ev$vectors))), e=5, label="eq3")

\[\begin{align} \label{eq3} A &= QDQ^{-1} \\ \begin{bmatrix} 1.000 & -3.000 \\ -2.000 & 4.000 \\ \end{bmatrix} &= \begin{bmatrix} 1.000 & -3.000 \\ -2.000 & 4.000 \\ \end{bmatrix} \end{align}\]

OK. Teraz przyda się funkcja korzystająca z wektorów i wartości własnych. Zauważmy

\[ A^n = \left( QDQ^{-1} \right)^n = QDQ^{-1}QDQ^{-1}\dots QDQ^{-1} \]

ponieważ \(Q^{-1}Q=I\) więc otrzymujemy

\[ A^n = \left( QDQ^{-1} \right)^n = QD^{n}Q^{-1} \]

co jest już o wiele prostsze w obliczeniach.

powm_e = function(M, n){
  ev = eigen(M)
  ev$vectors %*% (diag(ev$values))^n %*% inv(ev$vectors)
}

Metoda iloczynowa

Tylko w celach porównania stwórzmy funkcję iloczynową

powm_i = function(M, n){
  R = M
  for(i in 2:n){R = R %*% M}
  R
}

Sprawdzenie poprawności funkcji

Sprawdźmy czy nasze funkcje działają poprawnie.

n = 3
As = A %*% A %*% A 
Ad = powm_d(A, n)
Ae = powm_e(A, n)
Ai = powm_i(A, n)
xx("A_s = ", xm(As), "\\\\",
   "A_d = ", xm(Ad), "\\\\",
   "A_e = ", xm(Ae), "\\\\",
   "A_i = ", xm(Ai))

\[ A_s = \begin{bmatrix} 37.000 & -81.000 \\ -54.000 & 118.000 \\ \end{bmatrix} \\ A_d = \begin{bmatrix} 37.000 & -81.000 \\ -54.000 & 118.000 \\ \end{bmatrix} \\ A_e = \begin{bmatrix} 37.000 & -81.000 \\ -54.000 & 118.000 \\ \end{bmatrix} \\ A_i = \begin{bmatrix} 37.000 & -81.000 \\ -54.000 & 118.000 \\ \end{bmatrix} \]

Sprawdzenie wydajności

Macierz 2x2

Zacznijmy od macierzy 2x2.

set.seed(123)
k = 2
A = matrix(runif(k^2), k, k)

ev = eigen(A)
ev$values
## [1]  1.22641818 -0.05582325

Najpierw podnoszona do kwadratu

microbenchmark(powm_d(A, 2), powm_e(A, 2), powm_i(A, 2), times=100) %>% 
    autoplot()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

Teraz nieco większa potęga

microbenchmark(powm_d(A, 20), powm_e(A, 20), powm_i(A, 20), times=100) %>% 
    autoplot()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

I bardzo duża

microbenchmark(powm_d(A, 200), powm_e(A, 200), powm_i(A, 200), times=25) %>% 
    autoplot()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

I baaardzo duża

microbenchmark(powm_d(A, 2000), powm_e(A, 2000), powm_i(A, 2000), times=25) %>% 
    autoplot()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

Macierz 4x4

set.seed(123)
k = 4
A = matrix(runif(k^2), k, k)

ev = eigen(A)
ev$values
## [1]  2.3454389  0.7155885 -0.7052334 -0.1660017

Najpierw podnoszona do kwadratu

microbenchmark(powm_d(A, 2), powm_e(A, 2), powm_i(A, 2), times=100) %>% 
    autoplot()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

Teraz nieco większa potęga

microbenchmark(powm_d(A, 20), powm_e(A, 20), powm_i(A, 20), times=100) %>% 
    autoplot()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

I bardzo duża

microbenchmark(powm_d(A, 200), powm_e(A, 200), powm_i(A, 200), times=25) %>% 
    autoplot()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

I baaardzo duża

microbenchmark(powm_d(A, 2000), powm_e(A, 2000), powm_i(A, 2000), times=25) %>% 
    autoplot()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

Macierz 10x10

set.seed(123)
k = 10
A = matrix(runif(k^2), k, k)

ev = eigen(A)
ev$values
##  [1]  5.0651595+0.0000000i -0.4950848+0.5752976i -0.4950848-0.5752976i
##  [4]  0.6450438+0.3678797i  0.6450438-0.3678797i -0.7059716+0.0000000i
##  [7]  0.6138627+0.0000000i -0.5172094+0.0000000i  0.4269464+0.0000000i
## [10]  0.1739663+0.0000000i

Najpierw podnoszona do kwadratu

microbenchmark(powm_d(A, 2), powm_e(A, 2), powm_i(A, 2), times=100) %>% 
    autoplot()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

Teraz nieco większa potęga

microbenchmark(powm_d(A, 20), powm_e(A, 20), powm_i(A, 20), times=100) %>% 
    autoplot()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

I bardzo duża

microbenchmark(powm_d(A, 200), powm_e(A, 200), powm_i(A, 200), times=25) %>% 
    autoplot()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

I baaardzo duża

microbenchmark(powm_d(A, 2000), powm_e(A, 2000), powm_i(A, 2000), times=25) %>% 
    autoplot()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.