3 Matrix summary
Napravicemo funkciju ciji cilj je da odredi osnovna svojstva matrice, poput dimenzije, ranga, sopstvenih vrednosti, invertibilnosti, dijagonalnosti i slicno. Nakon toga napravicemo i funkciju koja stampa sazetak osnovnih svojstava matrice u lepom formatu, sa glue bibliotekom.
Pisacemo red po red funkcije, uz detaljna objasnjenja, pa na kraju dati kompletan kod.
3.1 Implementacija funkcije matrix_summary
Definisimo prvo funkciju, primace argument mat koji predstavlja matricu
Izlaz ce nam biti lista koja sadrzi trazena svojstva matrice. Biramo listu kao izlazni tip, jer zelimo da skladistimo objekte raznih dimenzija i tipova, pa se ne mozemo zadovoljiti vektorskim tipom. Dakle, kreiramo praznu listu, koju cemo dalje popunjavati.
Za pocetak dodacemo u listu izlaza matricu koju analiziramo, kao i dimenzije
Sledeci je rang matrice. Njega cemo racunati pomocu funkcije rankMatrix iz Matrix paketa
Da li je matrica kvadratna cemo proveriti uporedjivanjem dimenzija matrice. Ukoliko je razlika broja kolona i vrsta jednaka 0, to znaci da je matrica kvadratna.
Dalja svojstva matrice koja cemo gledati su: kvadratnost, inveritbilnost, inverz, sopstvene vrednosti, simetricnost i dijagonalnosti. Ova svojstva imaju smisla samo za kvadratne matrice, pa ce to biti uslov da razmatramo ova svojstva.
Da li je matrica singularna (tj. nema inverz) proveravamo tako sto pogledamo da li je rang matrice razlicit od broja vrsta, jer je matrica invertibilna akko je punog ranga.
Ukoliko matrica nije singularna, ima smisla izracunati joj inverz. Za to koristimo funkciju solve, koja sluzi za resavanje sistema oblika \(Ax = b\), dok ako joj se prosledi samo matrica, ona vraca inverz te matrice kao rezultat (za detalje pogledati dokumentaciju ?solve).
Za determinantu, imamo funkciju det…
Sopstvene vrednosti racunamo funkcijom eigen, koja kao rezultat vraca listu, ciji je jedan od elemenata $values – vektor sopstvenih vrednosti.
Simetricnost matrice proveravamo tako sto posmatramo da li je jednaka svom transponatu.
Za odredjivanje da li je matrica dijagonalna, implementiracemo u nastavku sopstvenu funkciju, za sad cemo samo dodati u listu svojstava njenu vrednost.
Na kraju zatvaramo if i vracamo konstruisanu izlaznu listu.
3.2 Provera dijagonalnosti
Radi provere dijagonalnosti, dacemo dve razlicite implementacije, pa cemo ih uporediti sa stanovista brzine i odabrti najbolju za nasu funkciju.
3.2.1 Implementacija 1
Prva implementacija se zasniva na jednostavnom prolazenju kroz matricu i ukoliko nadjemo element van dijagonale koji nije nula, vratimo FALSE.
is_diagonal <- function(mat) {
# Ukoliko argument mat nije matrica, javljamo gresku
if (!is.matrix(mat)) {
stop("Not a matrix")
}
# Prolazimo kroz sve redove...
for (i in 1:nrow(mat)) {
# ...i sve kolone
for (j in 1:ncol(mat)) {
# ako smo van dijagonale, i matrica na tom mestu nije nula,
# zakljucujemo da matrica nije dijagonalna, i vracamo FALSE
if (i != j && mat[i,j] != 0)
return(FALSE)
}
}
# U opstem slucaju vracamo TRUE
return(TRUE)
}3.2.2 Implementacija 2
Druga implementacija se oslanja na to da, kada dijagonalnu matricu element po element pomnozimo sa jedinicnom, na kraju opet dobijamo istu tu matricu. Par primera:
## [,1] [,2]
## [1,] 1 3
## [2,] 2 4## [,1] [,2]
## [1,] 3 0
## [2,] 0 8## [,1] [,2]
## [1,] 1 0
## [2,] 0 1## [,1] [,2]
## [1,] 1 0
## [2,] 0 4## [,1] [,2]
## [1,] 3 0
## [2,] 0 8To svojstvo mozemo lako implementirati…
3.2.3 Poredjenje brzine implementacija
Za poredjenje brzine implementacija koristicemo biblioteku microbenchmark.
Uporedicemo brzinu funkcija is_diagonal i is_diagonal2 za 2 razlicite matrice, jednu dijagonalnu i jednu ne-dijagonalnu.
diag_matrix <- diag(300) # 100x100 dijagonalna matrica
nondiag_matrix <- matrix(rnorm(9e4), ncol = 300) # random matrica 100x100Poredjenje za nedijagonalnu:
## Unit: microseconds
## expr min lq mean median uq
## is_diagonal(nondiag_matrix) 1.413 2.210 57.95029 3.7830 5.534
## is_diagonal2(nondiag_matrix) 363.019 371.921 611.81415 399.8815 761.445
## max neval
## 5395.846 100
## 3837.053 100U ovom slucaju je brza prva implementacija. Pogledajmo i poredjenje za dijagonalnu matricu.
## Unit: microseconds
## expr min lq mean median uq
## is_diagonal(diag_matrix) 6119.098 6202.453 6269.4888 6247.8910 6307.967
## is_diagonal2(diag_matrix) 545.611 911.607 948.7209 923.1235 945.001
## max neval
## 6595.381 100
## 3884.228 100U ovom slucaju je druga implementacija ubedljivo brza. Ako pogledamo pazljivije, druga implementacija ima u oba slucaja isto vreme, dok se u prvoj implementaciji sa for petljom vremena drasticno razlikuju.
Razlog za to je sto u is_diagonal, ako naidjemo na nenula element van dijagonale, odmah izlazimo iz funkcije, a to se desava prilicno brzo, jer vec u drugoj iteraciji petlje dodjemo na neki element van dijagonale. S druge strane, ako imamo dijagonalnu matricu, moramo da prodjemo kroz sve elemente matrice da bismo zakljucili da zaista nisu nule, a to zahteva veliko vreme.
Ali avaj! Kada razmislite, druga, brza implementacija uvek mnozi svaki element sa svakim, pa poredi da li su svi elementi neke dve matrice jednaki, sto je mnogo operacija.
Zasto je prolazak kroz sve elemente matrice u for petlji bio toliko skup, a u is_diagonal2, gde mnozimo sve elemente matrice (i cak imamo veci broj operacija!) nemamo tu kaznu u vremenu izvrsavanja?
Odgovor je u tome sto su mnozenje matrica, poredjenja matrica, i sve operacije nad vektorima i matricama u R-u u pozadini implmentirane u programskom jeziku C ili Fortran-u, koji su kompajlirani jezici i neuporedivo su brzi nego rucno prolazenje kroz petlju u R.
Dakle, glavno pravilo kod pisanja brzog koda u R-u je da se oslanjate u sto vecoj meri na ugradjenje funkcije, umesto pravljenja svojih. Takodje, for petlje treba izbegavati u sirokom luku, jer vrlo postoji bolje resenje, a petlje dovode do slabe brzine kode.
3.3 Funkcija print_matrix_summary
Sada cemo napisati jednostavnu funkciju koja ce nase zakljucke iz funkcije matrix_summary napisati kao kratak opis matrice. Neki primeri koje bismo ocekivali:
“m by n matrix of rank r” “n by n square, symmetric, nonsingular matrix of rank n” “n by n square, singular matrix of rank 7”
Definisimo prvo funkciju, primace argument mat koji predstavlja matricu
Ucitacemo biblioteku glue koju cemo koristiti za pravljenje poruke.
Ucitacemo rezultat matrix_summary funkcije, na osnovu kog cemo praviti sazeti opis matrice.
Imamo dva glavna slucaja - kvadratne i nekvadratne matrice.
Ako matrica nije kvadratna, nemamo neki poseban rezultat, vec samo znamo dimenziju i rang, pa cemo samo to istampati
# ako nije kvadratna...
if (!mat_summary$is_square) {
# ... istampaj...
print(
# ... poruku koja je oblika "m by n matrix of rank r"
# funkcija glue dopusta koriscenje vrednosti promenljivih ako se pisu
# unutar zagrada { } i on ce zameniti taj izraz njegovom vrednoscu.
glue("{dim(mat)[1]} by {dim(mat)[2]} matrix of rank {mat_summary$rank}")
)
}Ako je matrica kvadratna, onda imamo nekoliko svojstava koje zelimo da proverimo i da ispisemo. String koji cemo da ispisemo cemo praviti nadovezivanjem svakog svojstva jedno po jedno, pocevsi od praznog stringa.
else {
# string koji sadrzi svojstva matrice cemo inicijalizovati kao prazan
properties <- ""
# ako je matrica simetricna, dodacemo string ", symmetric" na string
# properties
if (mat_summary$is_symmetric)
properties <- glue(properties, ", symmetric")
# ako je matrica dijagonalna, dodajemo string ", diagonal" na string
# properties (kao rezultat dobijamo ", symmetric, diagonal")
if (mat_summary$is_diagonal)
properties <- glue(properties, ", diagonal")
# konacno, dodajemo string ", singular" ako je matrica singularna
# (nema inverz) ili string ", nonsingular" ako matrica nije singularna.
properties <- glue(properties,
ifelse(mat_summary$is_singular,
", singular",
", nonsingular"))
# properties string sada ima oblik slican ovom:
# ", symmetric, diagonal, nonsingular" i ideja je da ga nadovezemo na string
# "m by n square", i na rezultat dodamo " matrix of rank r", cime bi kao
# rezultat dobili string tipa:
# "m by n square"+", symmetric, diagonal, nonsingular"+" matrix of rank r" =
# "m by n square, symmetric, diagonal, nonsingular matrix of rank r"
print(
# funkcija glue po default-u nadovezuje argumente koji su joj prosledjeni.
glue(
"{dim(mat)[1]} by {dim(mat)[2]} square{properties}",
" matrix of rank {mat_summary$rank}"
)
)
}3.4 Celokupni kod
# Returns whether the given matrix is diagonal
is_diagonal <- function(mat) {
if (!is.matrix(mat)) {
stop("Not a matrix")
}
for (i in 1:nrow(mat)) {
for (j in 1:ncol(mat)) {
if (i != j && mat[i,j] != 0)
return(FALSE)
}
}
return(TRUE)
}
is_diagonal2 <- function(mat) {
if (!is.matrix(mat)) {
stop("Not a matrix")
}
all(mat * diag(nrow(mat)) == mat)
}
# Matrix summary
# Inputs: matrix mat
# Outputs: Basic properties od the matrix
# E.g. Square? Diagonal? Symmetric? Rank, dimensions,
# eigenvalues, inverse, Singular?, determinant
matrix_summary <- function(mat) {
output <- list()
output$matrix <- mat
output$dimensions <- dim(mat)
library(Matrix)
output$rank <- rankMatrix(mat)
output$is_square <- diff(dim(mat)) == 0
if (output$is_square) {
output$is_singular <- !(output$rank == dim(mat)[1])
if (!output$is_singular) {
output$inverse <- solve(mat)
}
output$determinant <- det(mat)
output$eigenvalues <- eigen(mat)$values
output$is_symmetric <- all(mat == t(mat))
output$is_diagonal <- is_diagonal2(mat)
}
class(output) <- "matrixSummary"
return(output)
}
# Prints out matrix description
# E.g. "m by n matrix of rank r"
# "n by n square, symmetric, nonsingular matrix of rank n"
# "n by n square, singular matrix of rank 7"
print_matrix_summary <- function(mat) {
library(glue)
mat_summary <- matrix_summary(mat)
if (!mat_summary$is_square) {
print(
glue("{dim(mat)[1]} by {dim(mat)[2]} matrix of rank {mat_summary$rank}")
)
} else {
properties <- ""
if (mat_summary$is_symmetric)
properties <- glue(properties, ", symmetric")
if (mat_summary$is_diagonal)
properties <- glue(properties, ", diagonal")
properties <- glue(properties,
ifelse(mat_summary$is_singular,
", singular",
", nonsingular"))
print(
glue(
"{dim(mat)[1]} by {dim(mat)[2]} square{properties}",
" matrix of rank {mat_summary$rank}"
)
)
}
}