Aplicació: Producte de matrius

En aquesta lliçò es veurà com calcular el producte de dues matrius i algunes de les seves aplicacions. A més, recordarem l’ús de la instrucció assert() que vam veure a la lliçò sobre el producte escalar.

Producte de matrius en matemàtiques

En matemàtiques, donades dues matrius $A$ i $B$, el primer que necessitem per definir el producte $AB$ és que la matriu $A$ tingui exactament el mateix nombre de columnes que el nombre de files de la matriu $B$. En aquest cas, si $A$ és una matriu $n \times m$ i $B$ una matriu $m \times r$, el seu producte serà una matriu $C$ de mida $n \times r$ tal que el seu element $c_{ij}$ (element que es troba a la fila $i$ i columna $j$) es calcula com:

$ c{ij} = \sum{k=1}^m a{ik}b{kj} $

Un exemple és el que podeu trobar al dibuix de dalt, on el producte d’una $2 \times 3$ per una $3 \times 2$ dóna com a resultat una matriu $2 \times 2$.

Producte de matrius en C++

Intentarem ara reproduir el procediment anterior amb llenguatge informàtic. Com ham dit, en primer lloc cal comprovar que les dimensions quadren: les columnes de A i les files de B han de coincidir. Això ho podem ssegurar amb la instrucció assert() que ja hem introduït. Així, la nostra funció podria comenár com es segueix, on suposem que la matriu A no serà buida:

using vectorD = vector<double>;
using matriuD = vector<vectorD>;

matriuD producte(const matriuD& A, const matriuD& B) 
{
    assert(A[0].size() == B.size()); // Atura el programa si les dimensions no coincideixen
    ...
}

Recordeu que, tot i que no modificarem les matrius A i B, les passarem per referència constant per millorar l’eficiència del programa. Pensem ara quina forma haurà de tenir el codi d’aquesta funció. Si anomenem les dimensions de les matrius com ho hem fet a l’apartat anterior, per recórrer la matriu resultat necessitarem dos bucles (de n i r iteracions, respectivament), i per calcular cadascun dels seus elements haurem de fer un altra bucle per calcular la suma dels m elements corresponents. Així, tindrem

using vectorD = vector<double>;
using matriuD = vector<vectorD>;

matriuD producte(const matriuD& A, const matriuD& B) 
{
    assert(A[0].size() == B.size()); // Atura el programa si les dimensions no coincideixen

    int n = A.size();
    int m = B.size();      // o tmabé m = A[0].size()
    int r = B[0].size();

    matriuD C(n, vectorD(r));               // Inicialitzem C a tot zeros
    for (int i = 0; i < n; ++i) {           // Recorre files de C i files de A
        for (int j = 0; j < r; ++j) {       // Recorre columnes de C i columnes de B
            for (int k = 0; k < m; ++k) {   // Recorre columnes de A i files de B
            C[i][j] += A[i][k] * B[k][j];
            }
        }
    }
    return C;
}

Potències de matrius quadrades

Suposem que tenim creada la funció producte que hem definit abans. Clarament, per elevar al quadrat una matriu no cal més que fer el producte per ella mateixa, com

matriuD A_quadrat = producte(A, A);

Fixeu-vos que en aquest cas l’assert que hi ha al principi de la funció equival a comprovar que la matriu A és quadrada.

Si el que volem ara és calcular potències més grans de la matriu A, ho podem fer aplicant una de les tècniques que ja hem estudiat: la recursivitat. Així, la següent funció recursiva calcularia la potència n-èssima de la matriu A, on suposarem que n és un nombre enter més gran o igual que $1$.

matriuD potencia(const matriuD& A, int n)
{
    if (n == 1) return A;                   // Cas base, A^1 = A
    return producte(A, potencia(A, n-1));   // A^n = A * A^(n-1)
}

Amb aquesta funció s’estan produïnt n productes de matrius. Com bé suposareu, fer un producte de matrius és una operació molt costosa (per la seva gran quantitat de bucles) així que ens interessa estalviar-nos quants més productes millor. És per això que en general farem servir un altre algorisme per calcular la potència n-èssima.

Fixem-nos en el següent: si $n$ és parell, la potència n-èssima es pot descomposar com el següent producte

$ A^n = A^{n/2} * A^{n/2} $

i si n és senar, no ens queda més remei que fer com abans i descomposar el producte com $A^n = A * A^{n-1}$. Però aleshores n-1 serà parell i podrem fer la descomposició d’abans! La gràcia d’aquesta descomposició és que, com que les dues matrius que multipliquem són iguals, no cal que la calculem dos cops, i així cada cop que ens trobem amb un n parell ens estarem estalviant la meitat dels productes que ens queden. Amb aquesta idea, la funció quedaria

matriuD potencia(const matriuD& A, int n)
{
    if (n == 1) return A;                   // Cas base, A^1 = A
    if (n % 2 == 0) {
        matriuD B = potencia(A, n/2);       // 😃 Només calculem B un cop
        return producte(B, B);
    }
    return producte(A, potencia(A, n-1));   // Si n és senar, ens conformem amb això
}

Així, el nombre de productes de matrius que estem fent ara serà de l’ordre de $\log(n)$.

Fibonacci amb potències de matrius

A la lliçó de recursivitat es van introduir els nombres de fibonacci i es va veure com implementar un algorisme recursiu per calcular-los. Recordem que els nombres de Fibonacci són $ 0,1,1,2,3,5,8,13,21,34,55,… $ on cada terme de la successió es calcula com la suma dels dos anteriors. Si denotem el n-èssim terme de la successió (començant per $0$) com $F_n$, la preopietat que compleixen aquests nombre es pot escriure com $Fn = F{n-1} + F_{n-2}$. Veurem ara com podem calcular $F_n$ gràcies a la tècnica antrior de les potències de matrius ràpides.

Considerem la matriu $ M = \left(\begin{array}{cc} 0 & 1 \ 1 & 1 \end{array}\right) = \left(\begin{array}{cc} F_0 & F_1 \ F_1 & F_2 \end{array}\right) $

Definirem també les matrius $M_n$ com $ M_n = \left(\begin{array}{cc} Fn & F{n+1} \ F{n+1} & F{n+2} \end{array}\right) $

Es pot veure clarament que $M = M_0$. Vegem que pasa si fem el producte de la matriu $M$ per la matriu $M_n$. Tenim que

$ \left(\begin{array}{cc} 0 & 1 \ 1 & 1 \end{array}\right) \left(\begin{array}{cc} F{n} & F{n+1} \ F{n+1} & F{n+2} \end{array}\right) = \left(\begin{array}{cc} F{n+1} & F{n+2} \ F{n}+F{n+1} & F{n+1}+F{n+2} \end{array}\right) = \left(\begin{array}{cc} F{n+1} & F{n+2} \ F{n+2} & F{n+3} \end{array}\right) $

És a dir, $M\cdot Mn = M{n+1}$ per qualsevol $n$! En particular tenim que $M^2 = M\cdot M_0 = M_1$, $M^3 = M\cdot M_1 = M2$, i així successivament. Per tant, $M^n = M{n-1}$ per qualsevol $n$. I ara ve la millor part: nosaltres sabem calcular la potència d’una matriu quadrada fent només uns $\log(n)$ productes 😃. Així, per trobar el terme $F_n$ de la successió de Fibonacci, només hem de calcular la potència $n+1$-èssima de $M$ (que és $M_n) i mirar el primer element, i ho podem fer en temps logarítmic (com que les matrius són $2\times 2$ podem considerar que el producte es fa en temps constant). Així, fent servir las funcions que hem anat creant al llarg de la lliçó (canviant tots els doubleper int en aquest cas), tindrem finalment:

int fibonacci(int n)
{
    matriuI M = {{0, 1}, {1, 1}};
    matriuI P = potencia(M, n+1);  // S'acomplirà sempre la hipòtesi n+1 >= 1.
    return P[0][0];
}

i podeu comprovar que d’aquesta manera el càlcul es fa molt més ràpid que amb la fórmula recursiva inicial. No obstant, els nombres de Fibonacci creixen ràpidament així que si voleu obtenir el seu valor per a valors de $n$ grans haureu de fer servir un tipus de dades capaç d’emagatzemar nombres més grans, com per exemple el long long.


Fòrum







Lliçons.jutge.org
Rafah Hajjar
Universitat Politècnica de Catalunya, 2019

Prohibit copiar. Tots els drets reservats.
No copy allowed. All rights reserved.