Inversione veloce di matrici: Strassen #2

In maniera molto simile a quanto è stato fatto per creare una funzione per eseguire velocemente moltiplicazioni tra grandi matrici con l'algoritmo di Strassen (vedi post precedente), scriviamo adesso funzioni per calcolare velocemente l'inversa di una matrice.

Per evitare di scrivere pagine e pagine di commenti e formule, come ho fatto per la moltiplicazione tra matrici, questa volta posto direttamente il codice delle funzioni create (il ragionamento che sta alla base è del tutto analogo). Copiate e incollate tutto il codice per poterlo visualizzare correttamente.

Funzione strassenInv(A)

strassenInv <- function(A){

 div4 <- function(A, r){
  A <- list(A)
  A11 <- A[[1]][1:(r/2),1:(r/2)]
  A12 <- A[[1]][1:(r/2),(r/2+1):r]
  A21 <- A[[1]][(r/2+1):r,1:(r/2)]
  A22 <- A[[1]][(r/2+1):r,(r/2+1):r]
  A <- list(X11=A11, X12=A12, X21=A21, X22=A22)
  return(A)
 }

        if (nrow(A) != ncol(A)) 
          { stop("only square matrices can be inverted") }

 is.wholenumber <-
     function(x, tol = .Machine$double.eps^0.5)  abs(x - round(x)) < tol

 if ( (is.wholenumber(log(nrow(A), 2)) != TRUE) || (is.wholenumber(log(ncol(A), 2)) != TRUE) )
   { stop("only square matrices of dimension 2^k * 2^k can be inverted with Strassen method") }

 A <- div4(A, dim(A)[1])

 R1 <- solve(A$X11)
 R2 <- A$X21 %*% R1
 R3 <- R1 %*% A$X12
 R4 <- A$X21 %*% R3
 R5 <- R4 - A$X22
 R6 <- solve(R5)
 C12 <- R3 %*% R6
 C21 <- R6 %*% R2
 R7 <- R3 %*% C21
 C11 <- R1 - R7
 C22 <- -R6
 
 C <- rbind(cbind(C11,C12), cbind(C21,C22))

 return(C)
}

Funzione strassenInv2(A)

strassenInv2 <- function(A){

 div4 <- function(A, r){
  A <- list(A)
  A11 <- A[[1]][1:(r/2),1:(r/2)]
  A12 <- A[[1]][1:(r/2),(r/2+1):r]
  A21 <- A[[1]][(r/2+1):r,1:(r/2)]
  A22 <- A[[1]][(r/2+1):r,(r/2+1):r]
  A <- list(X11=A11, X12=A12, X21=A21, X22=A22)
  return(A)
 }

 strassen <- function(A, B){
  A <- div4(A, dim(A)[1])
  B <- div4(B, dim(B)[1])
  M1 <- (A$X11+A$X22) %*% (B$X11+B$X22)
  M2 <- (A$X21+A$X22) %*% B$X11
  M3 <- A$X11 %*% (B$X12-B$X22)
  M4 <- A$X22 %*% (B$X21-B$X11)
  M5 <- (A$X11+A$X12) %*% B$X22
  M6 <- (A$X21-A$X11) %*% (B$X11+B$X12)
  M7 <- (A$X12-A$X22) %*% (B$X21+B$X22)

  C11 <- M1+M4-M5+M7
  C12 <- M3+M5
  C21 <- M2+M4
  C22 <- M1-M2+M3+M6
 
  C <- rbind(cbind(C11,C12), cbind(C21,C22))
  return(C)
 }

        if (nrow(A) != ncol(A)) 
          { stop("only square matrices can be inverted") }

 is.wholenumber <-
     function(x, tol = .Machine$double.eps^0.5)  abs(x - round(x)) < tol

 if ( (is.wholenumber(log(nrow(A), 2)) != TRUE) || (is.wholenumber(log(ncol(A), 2)) != TRUE) )
   { stop("only square matrices of dimension 2^k * 2^k can be inverted with Strassen method") }

 A <- div4(A, dim(A)[1])

 R1 <- strassenInv(A$X11)
 R2 <- strassen(A$X21 , R1)
 R3 <- strassen(R1 , A$X12)
 R4 <- strassen(A$X21 , R3)
 R5 <- R4 - A$X22
 R6 <- strassenInv(R5)
 C12 <- strassen(R3 , R6)
 C21 <- strassen(R6 , R2)
 R7 <- strassen(R3 , C21)
 C11 <- R1 - R7
 C22 <- -R6
 
 C <- rbind(cbind(C11,C12), cbind(C21,C22))

 return(C)
}

funzione strassenInv3(A)

strassenInv3 <- function(A){

 div4 <- function(A, r){
  A <- list(A)
  A11 <- A[[1]][1:(r/2),1:(r/2)]
  A12 <- A[[1]][1:(r/2),(r/2+1):r]
  A21 <- A[[1]][(r/2+1):r,1:(r/2)]
  A22 <- A[[1]][(r/2+1):r,(r/2+1):r]
  A <- list(X11=A11, X12=A12, X21=A21, X22=A22)
  return(A)
 }

 strassen <- function(A, B){
  A <- div4(A, dim(A)[1])
  B <- div4(B, dim(B)[1])
  M1 <- (A$X11+A$X22) %*% (B$X11+B$X22)
  M2 <- (A$X21+A$X22) %*% B$X11
  M3 <- A$X11 %*% (B$X12-B$X22)
  M4 <- A$X22 %*% (B$X21-B$X11)
  M5 <- (A$X11+A$X12) %*% B$X22
  M6 <- (A$X21-A$X11) %*% (B$X11+B$X12)
  M7 <- (A$X12-A$X22) %*% (B$X21+B$X22)

  C11 <- M1+M4-M5+M7
  C12 <- M3+M5
  C21 <- M2+M4
  C22 <- M1-M2+M3+M6
 
  C <- rbind(cbind(C11,C12), cbind(C21,C22))
  return(C)
 }

 strassen2 <- function(A, B){
  A <- div4(A, dim(A)[1])
  B <- div4(B, dim(B)[1])
  M1 <- strassen((A$X11+A$X22) , (B$X11+B$X22))
  M2 <- strassen((A$X21+A$X22) , B$X11)
  M3 <- strassen(A$X11 , (B$X12-B$X22))
  M4 <- strassen(A$X22 , (B$X21-B$X11))
  M5 <- strassen((A$X11+A$X12) , B$X22)
  M6 <- strassen((A$X21-A$X11) , (B$X11+B$X12))
  M7 <- strassen((A$X12-A$X22) , (B$X21+B$X22))

  C11 <- M1+M4-M5+M7
  C12 <- M3+M5
  C21 <- M2+M4
  C22 <- M1-M2+M3+M6

  C <- rbind(cbind(C11,C12), cbind(C21,C22))
  return(C)
 }

        if (nrow(A) != ncol(A)) 
          { stop("only square matrices can be inverted") }

 is.wholenumber <-
     function(x, tol = .Machine$double.eps^0.5)  abs(x - round(x)) < tol

 if ( (is.wholenumber(log(nrow(A), 2)) != TRUE) || (is.wholenumber(log(ncol(A), 2)) != TRUE) )
   { stop("only square matrices of dimension 2^k * 2^k can be inverted with Strassen method") }

 A <- div4(A, dim(A)[1])

 R1 <- strassenInv2(A$X11)
 R2 <- strassen2(A$X21 , R1)
 R3 <- strassen2(R1 , A$X12)
 R4 <- strassen2(A$X21 , R3)
 R5 <- R4 - A$X22
 R6 <- strassenInv2(R5)
 C12 <- strassen2(R3 , R6)
 C21 <- strassen2(R6 , R2)
 R7 <- strassen2(R3 , C21)
 C11 <- R1 - R7
 C22 <- -R6
 
 C <- rbind(cbind(C11,C12), cbind(C21,C22))

 return(C)
}

Proviamo a fare qualche prova. Innanzitutto controlliamo che la funzione esegua correttamente i calcoli, confrontandoli con il risultato della funzione standard di R (solve):

A <- matrix(trunc(rnorm(512*512)*100), 512,512)

all( round(solve(A),8) == round(strassenInv(A),8) )
[1] TRUE

all( round(solve(A),8) == round(strassenInv2(A),8) )
[1] TRUE

all( round(solve(A),6) == round(strassenInv3(A),6) )
[1] TRUE

La funzione esegue correttamente i calcoli. Esiste però un problema di approssimazione, per cui le prime due funzioni riportano il risultato identico fino alla ottava cifra decimale, mentre il secondo fino alla sesta. Più che un problema di calcoli, è un problema di espressione dei numeri in formato binario e 32 bit, che determina questi errori.

Ora analizziamo i tempi di esecuzione. Riporto in tabella il risultato ottenuto eseguendo il seguente codice:

Time computation

A <- matrix(trunc(rnorm(512*512)*100), 512,512)
system.time(solve(A))
system.time(strassenInv(A))
system.time(strassenInv2(A))
system.time(strassenInv3(A))

A <- matrix(trunc(rnorm(1024*1024)*100), 1024,1024)
system.time(solve(A))
system.time(strassenInv(A))
system.time(strassenInv2(A))
system.time(strassenInv3(A))

A <- matrix(trunc(rnorm(2048*2048)*100), 2048,2048)
system.time(solve(A))
system.time(strassenInv(A))
system.time(strassenInv2(A))
system.time(strassenInv3(A))

A <- matrix(trunc(rnorm(4096*4096)*100), 4096,4096)
system.time(solve(A))
system.time(strassenInv(A))
system.time(strassenInv2(A))
system.time(strassenInv3(A))

I risultati sono abbastanza evidenti, e utilizzando una modificazione dell'algoritmo di Strassen per l'inversione di matrici, esiste un effettivo risparmio di tempo.

Nuovamente valgono le due raccomandazioni già fatte:
- il codice è da migliorare; se qualcuno vuole aiutarmi, sarò lieto di aggiornare il mio codice
- se ritenete utili queste funzioni e le usate per qualche lavoro, una citazione è sempre gradita

Moltiplicazione veloce tra matrici: algoritmo di Strassen

Nel pdf qui sotto ho implementato una prima versione dell'algoritmo di Strassen in R, utile per svolgere rapidamente moltiplicazioni tra matrici.


Il codice è sicuramente da rivedere, e perfezionare, però può essere un utile inizio per una implementazione migliore.
Vorrei ringraziare l'utente G. Grothendieck per avermi suggerito un modo per creare matrici quadrate più grandi partendo da matrici più piccole (capirete leggendo il pdf), su StackOverFlow.

Se utilizzate il mio codice, è sempre cortese specificare la fonte. Inoltre se qualcuno vuole collaborare, perfezionando il codice, è assolutamente ben accetto!

Conversione numeri decimali in IEEE-754-32bit

Per un pò di teoria sullo standard IEEE-754 rimando alla pagina di Wikipedia. Qui posto solamente la funzione per effettuare la conversione in R.


Innanzitutto scriviamo alcune funzioni per eseguire la conversione dei numeri decimali in numeri binari:


decInt_to_8bit <- function(x, precs) {
q <- c()
r <- c()
xx <- c()
for(i in 1:precs){
xx[1] <- x
q[i] <- xx[i] %/% 2
r[i] <- xx[i] %% 2
xx[i+1] <- q[i]
}
rr <- rev(r)
return(rr)
}

devDec_to_8bit <- function(x, precs) {
nas <- c()
nbs <- c()
xxs <- c()
for(i in 1:precs)
{
xxs[1] <- x*2
nas[i] <- (xxs[i]) - floor(xxs[i])
nbs[i] <- trunc(xxs[i], 1)
xxs[i+1] <- nas[i]*2
}
return(nbs)
}


La prima funzione permette di convertire numeri interi (o prima della virgola) in codice binario, mentre la seconda permette di eseguire la conversione dei numeri decimali (dopo la virgola).
Il parametro precs indica il numero di bit da utilizzare.

Ad esempio:


decInt_to_8bit(11, 8)
[1] 0 0 0 0 1 0 1 1


Il numero 11, in codice binario a 8 bit è: 00001011.


devDec_to_8bit(0.625, 8)
[1] 1 0 1 0 0 0 0 0


Il numero 0.625 in codice binario a 8 bit è 10100000. Aumentando la precisione, nel primo caso permettiamo il calcolo di numeri superiori a 256 (limite massimo per gli 8-bit), mentre nel secondo caso permette di aumentare la precisione dopo la virgola:


devDec_to_8bit(0.3, 8)
[1] 0 1 0 0 1 1 0 0
devDec_to_8bit(0.3, 16)
[1] 0 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0


Ritornando ai primi due esempi, possiamo eliminare gli zeri in abbondanza con le seguenti due funzioni:


remove.zero.aft <- function(a) {
n <- length(a)
for(i in n:1){
if (a[n]==0) a <- a[-n]
else return(a)
n <- n-1
}
}

remove.zero.bef <- function(a) {
n <- length(a)
for(i in 1:n){
if (a[1]==0) a <- a[-1]
else return(a)
}
}


Abbiamo quindi:


remove.zero.bef(decInt_to_8bit(11, 8))
[1] 1 0 1 1

remove.zero.aft(devDec_to_8bit(0.625, 8))
[1] 1 0 1


Possiamo unire le due funzioni per crearne una che converte numeri decimali con la virgola in numeri binari:


dec.to.nbit <- function(x,n) {
aa <- abs(trunc(x, 1))
bb <- abs(x) - abs(trunc(x))

q <- c()
r <- c()
xx <- c()
for(i in 1:n){
xx[1] <- aa
q[i] <- xx[i] %/% 2
r[i] <- xx[i] %% 2
xx[i+1] <- q[i]
}
rr <- rev(r)

nas <- c()
nbs <- c()
xxs <- c()
for(i in 1:n)
{
xxs[1] <- bb*2
nas[i] <- (xxs[i]) - floor(xxs[i])
nbs[i] <- trunc(xxs[i], 1)
xxs[i+1] <- nas[i]*2
}

bef <- paste(remove.zero.bef(rr), collapse="")
aft <- paste(remove.zero.aft(nbs), collapse="")
bef.aft <- c(bef, aft)
strings <- paste(bef.aft, collapse=".")
return(strings)
}


Abbiamo quindi:


dec.to.nbit(11.625,8)
[1] "1011.101"


---

Imparato quindi a convertire i numeri decimali in numeri binari, possiamo passare alla conversione in IEEE-754 single precision (ossia 32bit), utilizzando la seguente funzione:


dec.to.ieee754 <- function(x) {
aa <- abs(trunc(x, 1))
bb <- abs(x) - abs(trunc(x))

rr <- decInt_to_8bit(aa, 32)

ppc <- 24 - length(remove.zero.bef(rr))

nbs <- devDec_to_8bit(bb, ppc)

bef <- remove.zero.bef(rr)
aft <- remove.zero.aft(nbs)

exp <- length(bef) - 1
mantissa <- c(bef[-1], aft)

exp.bin <- decInt_to_8bit(exp + 127, 16)
exp.bin <- remove.zero.bef(exp.bin)

first <- c()
if (sign(x)==1) first=c(0)
if (sign(x)==-1) first=c(1)

ieee754 <- c(first, exp.bin, mantissa, rep(0, 23-length(mantissa)))
ieee754 <- paste(ieee754, collapse="")

return(ieee754)
}


Il numero 11.625 in IEEE-754 è:


dec.to.ieee754(11.625)
[1] "01000001001110100000000000000000"


Per utilizzare la funzione dec.to.ieee754 servono le funzioni scritte sopra: decInt_to_8bit, devDec_to_8bit, remove.zero.bef, remove.zero.aft.

Per sicurezza potete confrontare questo output con quello che si ottiene da questo convertitore online

Regressione polinomiale #2

Ho già parlato della regressione polinomiale in un mio post precedente, dove ho focalizzato l'attenzione essenzialmente sui comandi da dare ad R, e sulla rappresentazione grafica dei modelli.
Ho deciso di scrivere un secondo post sullo stesso argomento per spiegare il procedimento formale per ottenere il modello, dando un forte rilievo ai calcoli manuali (più volte nel mio blog ho dato questa impostazione, per non ridurre la pratica statistica a semplici comandi testuali con R).


Nel pdf che ho scritto, vedremo quindi come calcolare formalmente i coefficienti beta del modello di regressione polinomiale. Per altro questa è una delle possibili applicazioni della decomposizione-LU di cui ho parlato nel mio ultimo post.

Potete accedere al file pdf, da questo link.

Decomposizione LU di una matrice

Il pane quotidiano della statistica è l'algebra lineare, e in particolare l'algebra delle matrici. Risulta pertanto fondamentale avere delle basi dei metodi utilizzati per poter apprendere a pieno i procedimenti.
In questo post voglio condividere un file pdf che ho scritto, per la soluzione di sistemi di equazioni con il metodo della decomposizione LU (o di Doolittle).

Ho cercato di descrivere tutti i passaggi nella maniera più chiara possibile, ma sicuramente una spiegazione "in diretta" risulta più utile. Pertanto ho selezionato i seguenti due video, disponibili su YouTube, che spiegano lo stesso procedimento che ho cercato di descrivere nel pdf.
Ho avuto la necessità di scrivere il pdf, in quanto nel video non viene spiegata come effettuare le permutazioni delle righe; procedimento che invece ho cercato di spiegare nel documento, in quanto R utilizza la decomposizione PLU (e non la decomposizione LU).

Da qui potete accedere al file LU-decompos.pdf.



Consiglio infine la lettura di questa pagina sullo stesso argomento.

NOTA: nel documento sarà presente qualche errore di sillabazione... la ragione è che non sono ancora riuscito a imparare come far sillabare correttamente le parole a latex!