Comprender un algoritmo para calcular un polinomio matricial

Question

Estoy tratando de entender esto algoritmo de Charles Van Loan para evaluar un polinomio matricial $p(\mathbf A)=\sum\limits_{k=0}^q b_k \mathbf A^k=b_0\mathbf I+b_1\mathbf A+\cdots$ (donde $\mathbf A$ est $n\times n$ ) :

$\displaystyle s=\lfloor q\rfloor, \quad r=\lfloor q/s\rfloor$
$\displaystyle \mathbf Y=\mathbf A^s$
$\displaystyle \text{for}\quad j=1,\dots,n$
$\displaystyle\quad\quad \mathbf y_0^{(j)}=\mathbf e_j$ ( $\mathbf e_j$ es el $j$ columna de la matriz identidad)
$\displaystyle\quad\quad \text{for}\quad k=1,\dots,s-1$
$\displaystyle\quad\quad\quad \mathbf y_k^{(j)}=\mathbf A\mathbf y_{k-1}^{(j)}$
$\displaystyle\quad\quad \mathbf f_0^{(j)}=\sum_{k=q}^{rs}b_k\mathbf y_{k-rs}^{(j)}$
$\displaystyle\quad\quad \text{for}\quad k=1,\dots,r$
$\displaystyle\quad\quad\quad \mathbf f_k^{(j)}=\mathbf Y\mathbf p_{k-1}^{(j)}+\sum_{i=0}^{s-1}b_{s(r-k)+i}\mathbf y_{i}^{(j)}$

He intentado repasarlo varias veces, pero parece que la cantidad (¿vector?) $\mathbf p_{k-1}^{(j)}$ no se definió en ninguna parte del documento. ¿Qué podría ser? Además, el documento afirma (o más probablemente es como yo (mal) entendí el documento) que sólo tres matrices deben ser almacenados para la evaluación: $\mathbf A$ , $\mathbf Y$ y el último $p(\mathbf A)$ . Sin embargo, no veo cómo organizar correctamente las cosas para que una matriz adicional para almacenar el $\mathbf y_k^{(j)}$ (de lo contrario, se necesitarían cuatro matrices en lugar de las tres que se indican en el artículo).

¿Podría alguien aclararme cómo aplicar correctamente este algoritmo? Mejor aún, ¿hay algún código disponible en algún sitio que implemente este algoritmo (soy agnóstico en cuanto al lenguaje, así que cualquier código que publique, probablemente pueda traducirlo al lenguaje que estoy usando actualmente).

Gracias por su ayuda.

Answer 1

El almacenamiento necesario para el $y^{(j)}_k$ es menor que una matriz "completa" en proporción como $s \ll n$ . Es decir, $s$ vectores de este tipo son necesarios para cada paso $j$ y el almacenamiento de una etapa puede reutilizarse en la siguiente.

La Pregunta afirma $s$ para ser el suelo de $q$ el grado del polinomio, pero esto es erróneo. Según la discusión del artículo que precede al Algoritmo I, $s$ puede ser cualquier número entero estrictamente comprendido entre 1 y $q$ . Aunque esto no garantiza $s \ll n$ lo permite. Véase la discusión de la aplicación del Algoritmo II al polinomio de Taylor para $e^A$ donde $s=4$ y $q=16$ fueron elegidos. [Obsérvese también la discusión sobre la minimización del "recuento de trabajo" en el Algoritmo II eligiendo $s$ en torno a $\sqrt{q/2}$ .]

Eso deja todavía el misterioso $p^{(j)}_k$ . Mi conjetura es que es un error tipográfico, que $f^{(j)}_k$ se calcula ajustando $Yf^{(j)}_{k-1}$ después de inicializar $f^{(j)}_0$ pero necesito hacer algunos garabatos para confirmar que funciona.

Añadido: Mi suposición anterior es correcta, que debemos leer $f^{(j)}_{k-1}$ en lugar de $p^{(j)}_{k-1}$ en un punto del Algoritmo II donde aparece este último.

El algoritmo puede entenderse como una combinación de dos ideas. Una es que el polinomio matricial $p(A)$ se evaluará columna a columna, es decir, aplicándola sucesivamente a los vectores unitarios $e_j$ , $j=1,..,n$ .

La otra idea es sencilla, aunque un poco engorrosa de explicar. Recordemos que $1 \lt s \lt q$ donde $q$ es el grado de $p$ . Sea el polinomio:

$$p(x) = \sum^q_{i=0} b_i x^i$$

se agruparán según la mayor potencia de $x^s$ que divide cada término y reescribir $p(x)$ en una expansión que implique potencias de $x^s$ y "coeficientes" que son polinomios de grado como máximo $s-1$ . Tomando $r = \lfloor q/s \rfloor$ y asumiendo los coeficientes indicados $b_m$ donde $m \gt q$ son cero:

$$p(x) = \sum^r_{k=0} x^{sk} \beta_k(x)$$

donde $\beta_k(x) = \sum^{s-1}_{i=0} b_{sk+i} x^i$ son los "coeficientes" del polinomio, como en el caso anterior.

Supongamos que $Y = A^s$ ya se ha calculado, por ejemplo, mediante exponenciación binaria u otro esquema.

A continuación, el bucle exterior del algoritmo II aplica iterativamente $p(A)$ a $n$ vectores unitarios $v = e_j$ como se ha explicado anteriormente. Para simplificar, nos limitaremos a describir la computación $p(A)v$ para un vector arbitrario $v$ de dimensión compatible, una columna de longitud $n$ donde $A$ est $n \times n$ . En efecto, omitimos el superíndice $(j)$ del resto del debate.

En primer lugar, un bucle para inicializar $y_0 = v$ y $y_i = A y_{i-1}$ para $i = 1,..,s-1$ .

El segundo es un bucle que esencialmente aplica el método de Horner para obtener $p(A)v$ utilizando la expansión en potencias de $A^s$ y "coeficientes" $\beta_k(A)v$ . Nótese que estos últimos son vectores columna, y para $k = 0,..,r$ evaluamos $\beta_k(A)v$ tomando la combinación lineal adecuada de los valores previamente inicializados $y_i$ . Es decir:

$$f_0 = \beta_r(A) v$$ $$f_k = Y f_{k-1} + \beta_{r-k}(A) v$$

para $k = 1,..,r$ . Esto demuestra que el término $p^{(j)}_{k-1}$ es un error tipográfico de $f^{(j)}_{k-1}$ en el Algoritmo II.

Answer 2

He aquí un código octave que implementa el "Algoritmo II" con bastante fidelidad. Las matrices nombradas con "_shifted" se indexan usando $i+1$ en el último componente siempre que el documento utilice $i,$ porque las matrices en matlab se indexan desde 1 en lugar de 0. Toda la aritmética se hace mod 29 para que sea fácil de comprobar contra una implementación de referencia más simple. No he tratado de hacerlo particularmente rápido; en realidad es sólo para demostrar que el algoritmo funciona.

q=20;
s=floor(sqrt(2*q));
r=floor(q/s);
n=100;
modulo=29;
A=randi(modulo,n,n);
b_shifted=randi(modulo,q+1,1);

pA_van_loan = zeros(n,n);

# Compute A^s
As=eye(n);
for k = 1:s
  As=A*As;
  As=mod(As, modulo);
end

# Van Loan algorithm II
for j = 1:n
  y_shifted = zeros(n,s);
  y_shifted(j,1) = 1;
  for k = 1:s-1
    y_shifted(:,k+1) = A*y_shifted(:,k-1+1);
    y_shifted(:,k+1) = mod(y_shifted(:,k+1), modulo);
  end

  f = zeros(n,1);
  for m = 0:(q-s*r)
    f = f + b_shifted(s*r+m+1) * y_shifted(:,m+1);
  end
  f = mod(f, modulo);
  for k = 1:r
    f = As * f;
    for m = 0:s-1
      f = f + b_shifted(s*(r-k)+m+1) * y_shifted(:,m+1);
    end
    f = mod(f, modulo);
  end

  pA_van_loan(:,j) = f;
end

# Compute polynomial naively for reference
pA_reference = zeros(n,n);
Aj = eye(n);
for j = 0:q
  pA_reference = pA_reference + b_shifted(j+1)*Aj;
  Aj=A*Aj;
  Aj=mod(Aj,modulo);
end
pA_reference=mod(pA_reference,modulo);

assert(pA_van_loan == pA_reference);

Answer 3

Aquí hay algo de código python que implementa el Algoritmo II de Van Loan. No está optimizado para la velocidad. En su lugar, he tratado de seguir los pasos del pseudocódigo en el papel tanto como sea posible. El código anterior puede ser optimizado y he mostrado en los comentarios dos posibilidades (pero estoy seguro de que hay más posibilidades).

import numpy as np

def van_loan(A: np.array, b: np.array) -> np.array:
    ''' Compute sum_{j=0}^q b[q] * A**q using Van Loan's algorithm II

    A needs to be a square matrix (there is no check in this function).
    b is a vector containing q+1 coefficients.
    '''
    n = A.shape[0]
    q = len(b)-1
    s = np.floor(np.sqrt(2*q)).astype(np.int)
    r = q // s
    Y = np.linalg.matrix_power(A, s)
    p_van_loan = np.zeros((n, n))
    y = np.zeros((n, s))
    for j in range(n):
        # Initialize y as being the j-th column of the n-by-n identity matrix
        y[:, 0] = 0
        y[j, 0] = 1

        # Compute the remaining y according to y_k = A*y_{k-1}
        for k in range(1, s):
            y[:, k] = np.dot(A, y[:, k-1])

        # Compute f_0
        f = np.zeros(n)
        for l in range(0, q-s*r+1):
            f += b[l+s*r]*y[:, l]
        # Compute f_0 alternative:
        # f = np.einsum("i,ni->n", b[s*r:q+1], y[:, :q-s*r+1])

        # Compute the remaining f iteratively
        for k in range(1, r+1):
            f = np.dot(Y, f)
            for l in range(0, s):
                f += b[s*(r-k)+l]*y[:, l]
            # This for loop can be replaced by the following (faster) 
            # line of code:
            # f += np.einsum("s,ns->n", b[s*(r-k):s*(r-k+1)], y[:, :s])

        # Set the j-th column to the computed f
        p_van_loan[:, j] = f

    return p_van_loan

Para demostrar que esto produce el mismo resultado que al aplicar la ecuación original de $p(\textbf{A})$ se puede utilizar el siguiente código:

def van_loan_reference(A, b):
    n = A.shape[0]
    Ap = np.eye(n)
    p = b[0]*Ap
    for j in range(1, len(b)):
        Ap = np.dot(Ap, A)
        p += b[j]*Ap
    return p

n = 3
q = 8
np.random.seed(0)
A = np.random.rand(n, n)
b = np.random.randn(q+1)
print(van_loan(A, b))
print(van_loan_reference(A, b))

Esto produce:

[[26.07842457 33.42897135 37.28842938]
 [22.01005978 31.30224965 33.21369951]
 [31.17386387 42.02787474 48.15436486]]
[[26.07842457 33.42897135 37.28842938]
 [22.01005978 31.30224965 33.21369951]
 [31.17386387 42.02787474 48.15436486]]