/
inv-ex2.Rmd
134 lines (99 loc) · 3.57 KB
/
inv-ex2.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
---
title: "Matrix inversion by elementary row operations"
author: "Michael Friendly"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{Matrix inversion by elementary row operations}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r nomessages, echo = FALSE}
knitr::opts_chunk$set(
warning = FALSE,
message = FALSE
)
options(digits=4)
par(mar=c(5,4,1,1)+.1)
```
The following examples illustrate the steps in finding the inverse of a matrix
using *elementary row operations* (EROs):
* Add a multiple of one row to another (`rowadd()`)
* Multiply one row by a constant (`rowmult()`)
* Interchange two rows (`rowswap()`)
These have the properties that they do not change the inverse. The method used here
is sometimes called the *Gauss-Jordan* method, a form of *Gaussian elimination*.
Another term is *(row-reduced) echelon form*.
Steps:
1. Adjoin the identity matrix to the right side of A, to give the matrix $[A | I]$
2. Apply row operations to this matrix until the left ($A$) side is reduced to $I$
3. The inverse matrix appears in the right ($I$) side
Why this works: The series of row operations transforms
$$ [A | I] \Rightarrow [A^{-1} A | A^{-1} I] = [I | A^{-1}]$$
If the matrix is does not have an inverse (is *singular*) a row of all zeros will appear
in the left ($A$) side.
### Load the `matlib` package
```{r }
library(matlib)
```
### Create a 3 x 3 matrix
```{r }
A <- matrix( c(1, 2, 3,
2, 3, 0,
0, 1,-2), nrow=3, byrow=TRUE)
```
### Join an identity matrix to A
```{r }
(AI <- cbind(A, diag(3)))
```
### Apply elementary row operations to reduce A to an identity matrix.
The right three cols will then contain inv(A).
We will do this three ways:
1. first, just using R arithmetic on the rows of `AI`
2. using the ERO functions in the `matlib` package
3. using the `echelon()` function
### 1. Using R arithmetic
```{r }
(AI[2,] <- AI[2,] - 2*AI[1,]) # row 2 <- row 2 - 2 * row 1
(AI[3,] <- AI[3,] + AI[2,]) # row 3 <- row 3 + row 2
(AI[2,] <- -1 * AI[2,]) # row 2 <- -1 * row 2
(AI[3,] <- -(1/8) * AI[3,]) # row 3 <- -.25 * row 3
```
Now, all elements below the diagonal are zero
```{r }
AI
#--continue, making above diagonal == 0
AI[2,] <- AI[2,] - 6 * AI[3,] # row 2 <- row 2 - 6 * row 3
AI[1,] <- AI[1,] - 3 * AI[3,] # row 1 <- row 1 - 3 * row 3
AI[1,] <- AI[1,] - 2 * AI[2,] # row 1 <- row 1 - 2 * row 2
AI
#-- last three cols are the inverse
(AInv <- AI[,-(1:3)])
#-- compare with inv()
inv(A)
```
### 2. Do the same, using matlib functions `rowadd()`, `rowmult()` and `rowswap()`
```{r }
AI <- cbind(A, diag(3))
AI <- rowadd(AI, 1, 2, -2) # row 2 <- row 2 - 2 * row 1
AI <- rowadd(AI, 2, 3, 1) # row 3 <- row 3 + row 2
AI <- rowmult(AI, 2, -1) # row 1 <- -1 * row 2
AI <- rowmult(AI, 3, -1/8) # row 3 <- -.25 * row 3
# show result so far
AI
#--continue, making above-diagonal == 0
AI <- rowadd(AI, 3, 2, -6) # row 2 <- row 2 - 6 * row 3
AI <- rowadd(AI, 2, 1, -2) # row 1 <- row 1 - 2 * row 2
AI <- rowadd(AI, 3, 1, -3) # row 1 <- row 1 - 3 * row 3
AI
```
### 3. Using `echelon()`
`echelon()` does all these steps *row by row*, and returns the result
```{r }
echelon( cbind(A, diag(3)))
```
It is more interesting to see the steps, using the argument `verbose=TRUE`. In
many cases, it is informative to see the numbers printed as fractions.
```{r }
echelon( cbind(A, diag(3)), verbose=TRUE, fractions=TRUE)
```