Multiplication et exponentiation de la chaîne matricielle

Si j'ai deux matrices et B , de dimensions 1000 × 2 et 2 × 1000 , respectivement, et que je veux calculer ( A B ) 5000 , il est plus efficace de réécrire d'abord l'expression sous la forme A ( B A ) 4999 B et ensuite seulement évaluer numériquement, car A B est de dimension 1000 × 1000 mais B A est de dimension 2 × 2 .$A$$B$$1000\times2$$2\times1000$$(AB)^{5000}$$A(BA)^{4999}B$$AB$$1000\times1000$$BA$$2\times2$

Je veux résoudre une version généralisée de ce problème. Existe-t-il un algorithme raisonnablement efficace (pas de force brute) pour optimiser une expression contenant:

• Variables matricielles libres de dimensions connues
• Produits de sous-expressions arbitraires
• Sous-expressions arbitraires élevées à la puissance naturelle

... de sorte qu'il faut le moins de travail pour évaluer numériquement, après avoir remplacé les variables matricielles libres par des valeurs matricielles concrètes?

Le problème de multiplication de la chaîne matricielle est un cas particulier de mon problème.

Éditer:

Ceci est une réponse provisoire. Cela me semble intuitivement juste, mais je n'ai aucune preuve que c'est correct. Si cela s'avère correct, je suis toujours intéressé par la preuve. (Si ce n'est pas correct, bien sûr, veuillez me corriger.)

Pour chaque produit élevé à une puissance, disons , considérons chaque permutation cyclique des facteurs:$(A_1 A_2 \ldots A_k)^n$

• $(A_1 A_2 \ldots A_k)^n$
• $A_1 (A_2 \ldots A_k A_1)^{n-1} A_2 \ldots A_k$
• $A_1 A_2 (A_3 \ldots A_k A_1 A_2)^{n-1} A_3 \ldots A_k$
• ...
• $A_1 A_2 \ldots A_{k-1} (A_k A_1 A_2 \ldots A_{k-1})^{n-1} A_k$

... récursivement. Chaque puissance doit être calculée en utilisant l'exponentiation par quadrature (évidemment), et tous les autres produits doivent être calculés en utilisant l'ordre optimal renvoyé par l'algorithme de multiplication de la chaîne matricielle.

Éditer:

L'idée décrite dans ma précédente édition n'est toujours pas optimale. L'algorithme d'exponentiation par quadrature évalue en fait des expressions de la forme ou A n K , où K n'est pas nécessairement la matrice d'identité. Mais mon algorithme ne considère pas la possibilité d'utiliser l'algorithme d'exponentiation par quadrature avec K différent de la matrice d'identité.$K A^n$$A^n K$$K$$K$

Avertissement: La méthode suivante n'a pas été rigoureusement prouvée comme étant optimale. Une preuve informelle est fournie.

Le problème se résume à trouver la commande la plus efficace lorsque l'on considère le carré du produit.

$(ABC)^{50}$$(ABC)^2$$ABCABC$$ABC$

$ABCABC$

$A(B(CA))BC$$A(B(CA))^{49}BC$

$(A_1 A_2 \cdots A_n)^m$$(A_1 A_2 \cdots A_n)^2$
$(A_1 A_2 \cdots A_n)^2$
$G$$A_1 \cdot A_2 \cdot G^{m-1} \cdot A_n$

$(AB)^n$$A$$B$$X \times Y$$Y \times X$$A$$B$

$X \times Y$
$Y \times X$
$Y \times Y$
$X \times X$

$X < Y$$Y ≤ X$

$X < Y$
$AB$$X \times X$$A$$B$$(AB)^n$

$Y ≤ X$
$BA$$Y \times Y$$A$$B$$A(BA)^{n-1}B$

$ABAB$

En utilisant plus de matrices, l'argument est similaire. Peut-être qu'une preuve inductive est possible? L'idée générale est que la résolution du MCM pour le carré trouvera la taille optimale pour les opérations avec toutes les matrices impliquées considérées.

Étude de cas:

julia> a=rand(1000,2);
julia> b=rand(2,1000);
julia> c=rand(1000,100);
julia> d=rand(100,1000);
julia> e=rand(1000,1000);

julia> @time (a*b*c*d*e)^30;
  0.395549 seconds (26 allocations: 77.058 MB, 1.58% gc time)

# Here I use an MCM solver to find out the optimal ordering for the square problem
julia> Using MatrixChainMultiply
julia> matrixchainmultiply("SOLVE_SQUARED", a,b,c,d,e,a,b,c,d,e)
Operation: SOLVE_SQUARED(A...) = begin  # none, line 1:
    A[1] * (((((A[2] * A[3]) * (A[4] * (A[5] * A[6]))) * (A[7] * A[8])) * A[9]) * A[10])
  end
Cost: 6800800

# Use the ordering found, note that exponentiation is applied to the group of 5 elements
julia> @time a*(((((b*c)*(d*(e*a)))^29*(b*c))*d)*e);
  0.009990 seconds (21 allocations: 7.684 MB)

# I also tried using the MCM for solving the problem directly
julia> @time matrixchainmultiply([30 instances of a,b,c,d,e]);
  0.094490 seconds (4.02 k allocations: 9.073 MB)

$A_1$$A_n$$A_i$$A_j$$O (n^3)$