5 Algebra Lineal

Muchos problemas reales en ciencia de datos, inteligencia artificial y aprendizaje automático implican el manejo y la manipulación de datos.

El álgebra lineal es fundamental en computación. El álgebra lineal proporciona herramientas para representar y operar con conjuntos de datos, como matrices y vectores, de manera eficiente.
Para técnicas de procesamiento de imágenes, se utilizan transformaciones lineales para realizar operaciones como rotaciones, escalados y traslaciones. Estas transformaciones se pueden llevar a cabo con álgebra lineal.
Muchos problemas en ciencia de la computación se pueden modelar mediante sistemas de ecuaciones lineales. La resolución de estos sistemas es esencial en diversas aplicaciones.
Los algoritmos de aprendizaje automático, especialmente en el ámbito de las redes neuronales y el deep learning, dependen en gran medida de operaciones matriciales.
En el diseño de algoritmos criptográficos, como RSA, el álgebra lineal es una herramienta fundamental.

En esta sección desarrollamos como trabajar con R para implementar:

las principales técnicas de álgebra numérica
los métodos básicos de resolución de sistemas de ecuaciones
los métodos de resolución de grandes sistemas de ecuaciones

Y en particular como trabajar con matrices dispersas que se hacen esenciales cuando manejamos cantidades ingentes de datos, para reducir la complejidad temporal y espacial de los algoritmos de extracción de conocimientos.

5.1 Algebra Lineal numérica en R

Resumimos en R las operaciones matriciales más importantes (ejecuta en tu ordenador para ver lo que hacen):

A <- matrix(c(1,1,0,2,3,1,0,1,0),nrow=3)
B <- matrix(c(-1,0,1,1,1,2,0,1,0),nrow=3)

A * B # multiplicación elemento a elemento
A %*% B  # Multiplicación normal de matrices
crossprod(A,B)  # A'B
crossprod(A)      # A'A
t(A)            # transpuesta de A
diag(A)         # extraer la diagonal - devuelve un vector
diag(1,4)       # devuelve la matriz identidad de orden 4
diag(c(1,3,1), 3) # devuelve la matriz diagonal con vector en la diagonal de orden 3
diag(diag(A)) # Ejecutar en R
diag(1, 3, 4)  # Ejecutar en R
solve(A) # Calcular la inversa
y<-eigen(A) # Calcular los autovalores y autovectores
y$val
y$vec

5.1.1 Paquete pracma

Es un paquete desarrollado para Algebra Lineal Numérica, optimización, ecuaciones diferenciales, series temporales, etc.

Instalar el paquete antes de usar los siguientes comandos (ejecuta en tu ordenador):

library(pracma)
# Usar estas matrices
A <- matrix(c(1,1,0,2,3,1,0,1,0),nrow=3)
v <- c(1,2,3) 
cond(A)  # devuelve el nº de condición de una matriz - debe estar cercano a 1
Polcar <- charpoly(A, info = TRUE) # polinomio característico de A
eye(3)    # devuelve la matriz identidad de orden 4
ones(3)    # devuelve la matriz de unos de orden 4
zeros(3)    # devuelve la matriz de ceros de orden 4
eye(3,4); ones(3,4); zeros(3,4)    # Ejecutar en R

Ver https://www.rdocumentation.org/packages/pracma/versions/1.9.9 para un listado de las funciones que integran pracma. Los usuarios de Matlab pueden ver muchas similitudes con los comandos usados en Matlab.

Es un paquete desarrollado para Algebra Lineal Numérica, optimización, ecuaciones diferenciales, series temporales, etc.

Instalar el paquete antes de usar los siguientes comandos (ejecuta en tu ordenador):

Diag(v, 1)    # devuelve la matriz diagonal con vector $v$  en la diagonal
Diag(A, 1)    # devuelve el vector de A que está en diagonal 1 
tril(A)   # devuelve la triangular inferior y la diagonal
triu(A)   # devuelve la triangular superior y la diagonal

hilb(3) # matriz de hilbert de orden 3
magic(3)  # matriz mágica de orden 3
Rank(A)  # el rango de una matriz
Trace(A) #la traza de una matriz

Para el cálculo de la inversa se usa el comando solve:

A <-matrix(rnorm(9),nrow=3); A1 <- solve(A)
Ident <- A%*%A1;Ident

             [,1]          [,2]        [,3]
[1,] 1.000000e+00  1.044384e-16 1.74335e-17
[2,] 1.929600e-16  1.000000e+00 1.44500e-17
[3,] 2.144307e-17 -4.914454e-17 1.00000e+00

# hay errores de redondeo - comprobación de igualdad exacta
 identical(A, solve(solve(A)))

[1] FALSE

# comprobación teniendo en cuenta la precisión de la máquina
 all.equal(A, solve(solve(A)))

[1] TRUE

5.1.2 Métodos Directos

Los métodos de resolución de SEL en R son los siguientes (Librería numérica LAPACK (http://www.netlib.org/lapack/) ).

solve
qr

A <-matrix(rnorm(9),nrow=3)
b <- c(1,2,-1)
x <- solve(A,b)
x
residuo <- A%*%x-b
norm(residuo,'1')
norm(residuo,'I')
norm(residuo,'2')

Dada una matriz \(A\) de tamaño \(m \times n\), el método QR factoriza la matriz como \(A=QR\) donde \(Q\) es una matriz ortogonal y \(R\) una matriz triangular superior. Se usa para resolver sistemas en los que \(m >n\) o los llamados sistemas sobredeterminados y/o mal condicionados.

A = matrix(c(1,   1,
             1,  -1,
            3,  -1,
            0, 1), nrow=4, byrow=TRUE)
b = c(2,0.1, 2,1)
b <- as.matrix(b)
Ab <- cbind(A,b) 
qr(A)$rank
qr(Ab)$rank
solve(A,b)# da error
x <- qr.solve(A,b)
x
residuo <- A%*%x-b
residuo
norm(residuo,'2')

5.2 Métodos Iterativos

Se usará el paquete pracma. Este paquete tiene una función que permite resolver métodos iterativos por los métodos de Jacobi, Gauss-Seidel y Richardson.

itersolve <- function(A, b, x0 = NULL,nmax = 1000, 
                tol = .Machine$double.eps^(0.5),
                method = c("Gauss-Seidel", "Jacobi",
                         "Richardson"))
....

Es más compleja aunque más potente en cuanto a control de errores de los parámetros de entrada, etc.

5.2.1 Jacobi

Para resolver un sistema por Jacobi (el procedimiento sería igual en Gauss-Seidel), ejecutaríamos lo siguiente:

1000 iteraciones
partimos del vector inicial nulo

library(pracma)
A=matrix(c(2, 3, 4, 1 ,-3, 0,1, -1, 4),nrow=3,byrow=TRUE)
b=c(-2, 4,6)
sol <-itersolve(A,b,c(0,0,0),nmax=1000,method=c("Jacobi"))
sol

Calcular el vector residuo (el error cometido)
Calcular la norma del vector residuo

5.3 Matrices Dispersas

Se llama matriz dispersa a matrices de gran tamaño con la mayoría de elementos nulos. Las más interesante matrices dispersas son las que se llaman estructuradas: Los elementos no nulos forman patrones regulares que son fáciles de introducir en la matriz.

Las estructuras de datos se optimizan para almacenar estos datos de forma mucho más eficiente y reducir el tiempo de acceso a los elementos.

Aplicaciones:

Análisis de Redes Sociales, minería de textos - procesamiento del lenguaje natural. Representando en una matriz por filas los documentos y por columnas las palabras.
Sistemas de recomendación. Representando en una matriz por filas los usuarios y por columnas las películas o canciones vistas o escuchadas por los usuarios.
Market basket analysis: usuarios versus productos.
Dinámica de fluidos - ecuaciones diferenciales
Criptografía
Mecánica estructural
Google page rank
Visión ordenador
En general en problemas en los que el tamaño es crítico: bigdata

5.3.1 Almacenamiento

El uso de matrices dispersas consigue una optimización del tamaño usado para representar la matriz. Ejemplo de la diferencia de cantidad de memoria necesitada:

library('Matrix')
M1 <- Matrix(0, nrow = 1000, 
                  ncol = 1000, 
                  sparse = TRUE)
M1[1,100]<-1; M1[2,100]<-1
M2 <- matrix(0, nrow = 1000, 
                  ncol = 1000)
M2[1,100]<-1; M2[2,100]<-1
print (object.size(M1))

5816 bytes

print (object.size(M2))

8000216 bytes

La matrix dispersa solo almacena los valores no nulos. Hay distintas formas de almacenar la información no nula.

El método más sencillo es por columnas. Los elementos se almacenan mediante tuplas con tres variables: \((i,j,A_{i,j})\) donde \(i,j\) indican la fila y la columna y \(A_{i,j}\) el valor no nulo.

M1 <- Matrix(nrow = 3, ncol = 3, data = 0, sparse = TRUE)
M <- as(M1,"TsparseMatrix") 
M[1,1] <- 1
M[2,2] <- 1
M[3,1] <- 1
M

3 x 3 sparse Matrix of class "dgTMatrix"
          
[1,] 1 . .
[2,] . 1 .
[3,] 1 . .

str(M)

Formal class 'dgTMatrix' [package "Matrix"] with 6 slots
  ..@ i       : int [1:3] 0 1 2
  ..@ j       : int [1:3] 0 1 0
  ..@ Dim     : int [1:2] 3 3
  ..@ Dimnames:List of 2
  .. ..$ : NULL
  .. ..$ : NULL
  ..@ x       : num [1:3] 1 1 1
  ..@ factors : list()

filas <- M@i+1
filas

[1] 1 2 3

columnas <- M@j+1
columnas

[1] 1 2 1

valores <- M@x
valores

[1] 1 1 1

El siguiente método sería el llamado Esquema CSR (Compressed sparse row): dos vectores para cada fila - un vector \(J\) que contiene los índices de esa fila y otro \(V\), que contiene los valores de la fila.

Usaremos los métodos implementados en el paquete Matrix:

v <- sample(c(0,1),25,
            replace = TRUE)
v

 [1] 0 0 0 0 1 0 0 1 0 0 1 0 0 0 1 0 1 0 0 0 0 0 1 0 0

M1 <- matrix(v, nrow = 5, 
                  ncol = 5)
M1

     [,1] [,2] [,3] [,4] [,5]
[1,]    0    0    1    0    0
[2,]    0    0    0    1    0
[3,]    0    1    0    0    1
[4,]    0    0    0    0    0
[5,]    1    0    1    0    0

M1_sparse <- as(M1, 
                "sparseMatrix")
M1_sparse

5 x 5 sparse Matrix of class "dgCMatrix"
              
[1,] . . 1 . .
[2,] . . . 1 .
[3,] . 1 . . 1
[4,] . . . . .
[5,] 1 . 1 . .

R permite almacenar y manipular de forma eficiente las matrices dispersas. Las operaciones habituales de matrices respetan la estructura dispersa:

M1_sparse*2

5 x 5 sparse Matrix of class "dgCMatrix"
              
[1,] . . 2 . .
[2,] . . . 2 .
[3,] . 2 . . 2
[4,] . . . . .
[5,] 2 . 2 . .

M1_sparse%*%M1_sparse

5 x 5 sparse Matrix of class "dgCMatrix"
              
[1,] . 1 . . 1
[2,] . . . . .
[3,] 1 . 1 1 .
[4,] . . . . .
[5,] . 1 1 . 1

El paquete Matrix contiene unas 70 clases de matrices distintas y miles de métodos para trabajar con ellas. Se usa el comando spMatrix para introducir matrices dispersas en memoria.

 A <- spMatrix(10,20, i = c(1,3:8), j = c(2,9,6:10),  x = 7 * (1:7)); A
str(A)
A[2:7, 12:20] <- rep(c(0,0,0,(3:1)*30,0), length = 6*9); A
# Si en algún cálculo deja de ser sparse, lo hace automáticamente
A + 1
3*A
# Operaciones matriciales con matrices dispersas
t(A) #matriz traspuesta
crossprod(A)  # t(A) %*% A
tcrossprod(A)  # A %*% t(A)
# Obtener una imagen de matriz dispersa
 image(A)

Algunos ejemplos más de cómo trabajar con matrices dispersas comparando con las mismas matrices pero densas.

 library('Matrix')
M1 <- matrix(0, nrow = 1000, ncol = 1000)
M2 <- Matrix(0, nrow = 1000, ncol = 1000, sparse = TRUE)
object.size(M1); object.size(M2)
#M1 no cambiará de tamaño si añadimos un elemento.
M1[500, 500] <- 1; M2[500, 500] <- 1
object.size(M1); object.size(M2) 
# Operaciones con matrices dispersas son realizadas muy eficientemente
M2 %*% rnorm(1000)
M2 + 2*M2
M2 - M2
t(M2) 
#Pueden usarse cBind, rBind
M3 <- cBind(M2, M2)
nrow(M3); ncol(M3)
M4 <- rBind(M2, M2)
nrow(M4); ncol(M4)

Convertir a sparse:

 library('Matrix')
M <- Matrix(1:30, 5, 6)
M[, c(1,3:5)] <- 0
M[c(1,3), ] <- 0
M
sM <- as(M, "sparseMatrix")
sM
# operaciones que se pueden hacer
3*sM
1+sM
sM>5

5.3.2 Matrices triangulares

Muchos problemas prácticos (resolución de sistemas de ecuaciones diferenciales por ejemplo) involucran resolución de sistemas de ecuaciones que tienen patrones particulares. El más sencillo es cuando el sistema se puede representar con matrices triangulares.

Construimos los tres vectores que forman la diagonal, se forma la matriz correspondiente y se resuelve con el comando solve.

library(pracma)
n <-5000 # Size of matrix
M <- matrix(n,n,data=0) # Creamos matriz cuadrada de ceros
# Unos en la diagonal
dinf <- rep(1,n-1)
dinf <- dinf[-1] # Trim aa to fit below the diagonal
dinf[n-1]  <-  0 # Excepto la primera y la última

# Creamos la diagonal principal
dprin <- rep(-1.99,n)
dprin[1] <- dprin[n] <- 1
# Creamos la diagonal superior
dsup <- rep(1,n-1)
dsup[n-1]  <- 0 # Except first and last element

Construimos los tres vectores que forman la diagonal, se forma la matriz correspondiente y se resuelve con el comando solve.

# Vector de los términos independientes
b  <-  rep(0,n)
b[1]  <-  0.5
b[n]  <-  0.69

# Assemble matrix
M[1,1:2]  <-  c(dprin[1],dsup[1])
M[n,(n-1):n]  <-  c(dinf[n-1],dprin[n])
for (i in 2:(n-1)) M[i,(i-1):(i+1)]  <-  c(dinf[i-1],dprin[i],dsup[i])
options(digits=3)
# Solve
system.time(soln  <-  solve(M,b))
#soln
res <-norm(M%*%soln-b,'1')
res

La solución anterior necesita construir la matriz completa. El paquete limSolve tiene varias funciones para resolver eficientemente sistemas de este tipo. No se almacenan las diagonales en una matriz con el consiguiente ahorro en espacio y en tiempo de computación.

# Solve.tridiag from the limSolve package. 
# This algorithm goes as O(n)- solve O(n^3).
require(limSolve)

n <-5000 # Size of matrix
M <- matrix(n,n,data=0) # Creamos matriz cuadrada de ceros
# Unos en la diagonal
dinf <- rep(1,n-1)
dinf <- dinf[-1] # Trim aa to fit below the diagonal
dinf[n-1]  <-  0 # Excepto la primera y la última

# Creamos la diagonal principal
dprin <- rep(-1.99,n)
dprin[1] <- dprin[n] <- 1

# Creamos la diagonal superior
dsup <- rep(1,n-1)
dsup[n-1]  <- 0 # Except first and last element

# Vector de los términos independientes
b  <-  rep(0,n)
b[1]  <-  0.5
b[n]  <-  0.69

system.time(tri.soln <- Solve.tridiag(dsup,dprin,dinf,b))

5.3.3 Matrices con bandas - Banded Matrix

El paquete limSolve permite resolver sistemas cuyas matrices tienen bloques o patrones determinados. Por ejemplo, permite resolver sistemas cuya matriz de los coeficientes es una matriz de bandas (bandex matrix).

require(limSolve)
options(digits=3)
set.seed(333)

n <- 5000  
# Lower diagonals
dn1 <- runif(n-1)
dn2 <- runif(n-2)

# Diagonal
bb <- runif(n)
# Upper diagonals
up1 <- runif(n-1)
up2 <- runif(n-2)

# Assemble matrix
abd <- rbind(c(0,0,up2),c(0,up1),bb,c(dn1,0),c(dn2,0,0))

B <- runif(n) # rhs of system

system.time(Band <- Solve.banded(abd, nup=2, nlow=2,B))

5.3.4 Matrices por bloques - Block Matrix

Es una matriz donde los elementos no nulos se encuentran en bloques de distintos tamaños y normalmente cerca de la diagonal.

# Define matrix dimensions, set elements to 0
AA = matrix (nr= 12, nc=12, 0)
# Enter matrix elements
AA[1,1:4] = c( 0.0, -0.98, -0.79, -0.15)
AA[2,1:4] = c(-1.00, 0.25, -0.87, 0.35)
AA[3,1:8] = c( 0.78, 0.31, -0.85, 0.89, -0.69, -0.98, -0.76, -0.82)
AA[4,1:8] = c( 0.12, -0.01, 0.75, 0.32, -1.00, -0.53, -0.83, -0.98)
AA[5,1:8] = c(-0.58, 0.04, 0.87, 0.38, -1.00, -0.21, -0.93, -0.84)
AA[6,1:8] = c(-0.21, -0.91, -0.09, -0.62, -1.99, -1.12, -1.21, 0.07)
AA[7,5:12] = c( 0.78, -0.93, -0.76, 0.48, -0.87, -0.14, -1.00, -0.59)
AA[8,5:12] = c(-0.99, 0.21, -0.73, -0.48, -0.93, -0.91, 0.10, -0.89)
AA[9,5:12] = c(-0.68, -0.09, -0.58, -0.21, 0.85, -0.39, 0.79, -0.71)
AA[10,5:12] = c( 0.39, -0.99, -0.12, -0.75, -0.68, -0.99, 0.50, -0.88)
AA[11,9:12] = c( 0.71, -0.64, 0.0, 0.48)
AA[12,9:12] = c( 0.08, 100.0, 50.00, 15.00)

 AA

Es una matriz donde los elementos no nulos se encuentran en bloques de distintos tamaños y normalmente cerca de la diagonal.

# Define matrix dimensions, set elements to 0
B = c(-1.92,-1.27,-2.12,-2.16,-2.27,-6.08,-3.03,-4.62,-1.02,
-3.52,0.55,165.08)

# matrix AA is divided into blocks as follows:
Top = matrix(nr=2, nc=4, data=AA[1:2,1:4])
Top
Bot = matrix(nr=2, nc=4, data=AA[11:12,9:12])
Bot
Blk1 = matrix(nr=4, nc=8, data=AA[3:6,1:8])
Blk1
Blk2 = matrix(nr=4, nc=8, data=AA[7:10,5:12])
Blk2

AR = array(dim=c(4,8,2),data=c(Blk1,Blk2))
# overlap is the sum of the number of rows of Top and Bot.
Solve.block(Top,AR,Bot,B,overlap=4)

5.3.5 Paquete arules

En muchas aplicaciones aparecen matrices dispersas. El paquete arules usa las matrices dispersas como base de su estructura de datos.

Definir una función denominada extract.candidates tal que reciba como entrada un conjunto de implicaciones obtenidas con apriori estableciendo confianza 1. Explicación de lo que hará la rutina:
- Una regla i es candidata si la parte izquierda de i está incluida en la parte izquierda de alguna regla j (\(i\neq j\)).
- La rutina devuelve el ordinal de toda las reglas candidatas del conjunto original.

Ayuda:

Usaré las siguientes funciones de arules: is.subset, lhs.
Coger un dataset el que querías y colocar primero soporte adecuado para tener pocas reglas.

library(arules)


Attaching package: 'arules'

The following objects are masked from 'package:base':

    abbreviate, write

data(Adult)
r1 <- apriori(Adult,parameter = list(confidence=1,support=0.135))

Apriori

Parameter specification:
 confidence minval smax arem  aval originalSupport maxtime support minlen
          1    0.1    1 none FALSE            TRUE       5   0.135      1
 maxlen target  ext
     10  rules TRUE

Algorithmic control:
 filter tree heap memopt load sort verbose
    0.1 TRUE TRUE  FALSE TRUE    2    TRUE

Absolute minimum support count: 6593 

set item appearances ...[0 item(s)] done [0.00s].
set transactions ...[115 item(s), 48842 transaction(s)] done [0.02s].
sorting and recoding items ... [23 item(s)] done [0.00s].
creating transaction tree ... done [0.01s].
checking subsets of size 1 2 3 4 5 6 7 8 done [0.01s].
writing ... [12 rule(s)] done [0.00s].
creating S4 object  ... done [0.00s].

extract.candidates <- function(reglas){
    matriz <- is.subset(lhs(reglas),lhs(reglas))
    colnames(matriz) <- 1:length(r1)
    rownames(matriz) <- 1:length(r1)
    diag(matriz) <- 0
    filas <- as.logical(rowSums(matriz))
    num.filas <- colnames(matriz)[filas]
return(num.filas)
}#end function

cand <- extract.candidates(r1)
cand

[1] "1" "2" "3" "4" "5" "6" "7" "8"

Las reglas son matrices dispersas:

izquierda <- lhs(r1[1])
izquierda

itemMatrix in sparse format with
 1 rows (elements/transactions) and
 115 columns (items)

inspect(izquierda)

    items                                           
[1] {relationship=Husband, hours-per-week=Over-time}

class(izquierda@data)

[1] "ngCMatrix"
attr(,"package")
[1] "Matrix"

#ejecuta en consola 
# izquierda@data

Vemos como ácceder a objetos que implementan la estructura de matriz dispersa.

library(Matrix)
v <- 1 + 1:15
M <- Matrix(v, 60, 50)
# ejecuta M en consola
str(M)

Formal class 'dgeMatrix' [package "Matrix"] with 4 slots
  ..@ Dim     : int [1:2] 60 50
  ..@ Dimnames:List of 2
  .. ..$ : NULL
  .. ..$ : NULL
  ..@ x       : num [1:3000] 2 3 4 5 6 7 8 9 10 11 ...
  ..@ factors : list()

# Ejemplo de operaciones con matrices dispersas
M1 <- M %*% t(M)
# Añadir filas o columnas
M2 <- cbind(1, M1)
#M2
M3 <- rbind(1,M2)
#M3
dim(M2)

[1] 60 61

dim(M3)

[1] 61 61

# Transformamos la matriz M3 añadiendo muchos 0's
M3[,3:55] <- 0
M3[1:15,2] <-0
# M3
#Comprobamos tamaños
object.size(M3)

30944 bytes

M4 <-as(M3, "sparseMatrix")
object.size(M4)

7424 bytes

# Comprobar si es dispersas
is(M3,"sparseMatrix")

[1] FALSE

is(M4,"sparseMatrix")

[1] TRUE

#Comprobamos si son identicas
identical(M3,M4)

[1] FALSE

#Resumen de ambas
summary(M3)

   Length     Class      Mode 
     3721 dgeMatrix        S4

#summary(M4)

# Operaciones aritméticas - deja de ser dispersa
M5 <- M4+1
is(M5,"sparseMatrix")

[1] FALSE

object.size(M5)

30944 bytes

# Operadores lógicos con matrices dispersas
# Mindex es una matriz dispersa lógica
Mindex <- (M4>=1)
# Mindex
Mindex[1,1]

[1] TRUE

Mindex[1,2]

[1] FALSE

# summary(Mindex)

Para acceder a los slots que proporciona la clase se usa ***

M10 <- Matrix(1,3,3,sparse = TRUE)
str(M10)

Formal class 'dsCMatrix' [package "Matrix"] with 7 slots
  ..@ i       : int [1:6] 0 0 1 0 1 2
  ..@ p       : int [1:4] 0 1 3 6
  ..@ Dim     : int [1:2] 3 3
  ..@ Dimnames:List of 2
  .. ..$ : NULL
  .. ..$ : NULL
  ..@ x       : num [1:6] 1 1 1 1 1 1
  ..@ uplo    : chr "U"
  ..@ factors : list()

# Dimensión de matríz
M10@Dim

[1] 3 3

# the non-zero elements of the matrix
M10@x

[1] 1 1 1 1 1 1

#uplo:
# upper triangle ("U") 
# lower triangle ("L") 
M10@uplo

[1] "U"

Al operar con matrices dispersas podrían aparecer 0’s en dicha matriz, en ese caso usar comando drop0.

# los valores de Mindex antes de hacer drop0
str(Mindex@x)

 logi [1:473] TRUE TRUE TRUE TRUE TRUE TRUE ...

Mindex <- drop0(Mindex)
# los valores de Mindex después
str(Mindex@x)

 logi [1:473] TRUE TRUE TRUE TRUE TRUE TRUE ...

Algunas operaciones con estas matrices:

sym <- is(M5, "symmetricMatrix")
sym

[1] FALSE

# Operaciones habituales trabajan con matrices dispersas
#exp(M5)
#sin(M5)
#round(M5,digits = 2)
#min(M5)
#sum(M5)
# Extraer índices y valores que cumplen ciertas condiciones
which(M5>10000)

 [1] 3493 3508 3523 3538 3553 3554 3568 3569 3583 3584 3598 3599 3613 3614 3615
[16] 3628 3629 3630 3643 3644 3645 3658 3659 3660 3673 3674 3675 3676 3688 3689
[31] 3690 3691 3703 3704 3705 3706 3718 3719 3720 3721

M5[M5>10000]

 [1] 10401 10401 10401 10401 10501 11201 10501 11201 10501 11201 10501 11201
[13] 10501 11251 12001 10501 11251 12001 10501 11251 12001 10501 11251 12001
[25] 10401 11201 12001 12801 10401 11201 12001 12801 10401 11201 12001 12801
[37] 10401 11201 12001 12801

5.3.6 Visualización

# Albers' example Matrix with "Difficult" Eigen Factorization
data(CAex)
dim(CAex)

[1] 72 72

print(image(CAex, main = "image(CAex)"))

5.3.7 SparseM

https://www.jstatsoft.org/article/view/v008i06

SparseM proporciona funciones básicas de R para álgebra lineal con matrices dispersas. El uso del paquete se ilustra mediante una familia de funciones de ajuste de modelos lineales que implementan métodos de mínimos cuadrados para problemas con matrices de diseño dispersas. Las aplicaciones que utilizan matrices dispersas de gran tamaño permiten mejorar considerablemente el rendimiento en cuanto a utilización de memoria y velocidad de cálculo. Más especificamente: estimación de parámetros en modelos lineales muy grandes; tablas de contingencia dispersas de gran tamaño

5.4 Ejercicios

5.4.1 1. Métodos Directos

Resolver el sistema \(Ax=b\) usando R.

5.4.2 2. Métodos Iterativos

istema \(Ax=b\) usando R.

Dadas en R la siguiente matriz \(A\), y los vectores \(b\) y \(x_0\), resolver \(Ax=b\) por Jacobi y Gauss-Seidel. Comprobar qué método funciona mejor (error, convergencia).

A = matrix(c(4, -1, -1, 0, 0, 0,
-1, 4, 0, -1, 0, 0,
-1, 0, 4, -1, -1, 0,
0, -1, -1, 4, 0, -1,
0, 0, -1, 0, 4, -1,
0, 0, 0, -1, -1, 4), nrow=6, byrow=TRUE)

b = c(1, 5, 0, 3,1,5)

x0 = c(0.25, 1.25, 0, 0.75, 0.25, 1.25)

Calcular el vector residuo (el error cometido)
Calcular la norma del vector residuo

5.4.3 3. Función con matriz dispersa

Definir una función tal que tome como entrada una matriz dispersa y devuelva:

\[degree.sparsity(M)=1-\frac{count.values.non.zero(A)}{total.posiciones(A)}\]

5.4.4 4. Sistema triangular

Resolver el siguiente sistema triangular \(Ax=b\) (dimensión n=10000),