2016年11月24日星期四

Use R to solve linear euqations

https://cran.r-project.org/web/packages/matlib/vignettes/linear-equations.html

Michael Friendly

2016-09-16

This vignette illustrates the ideas behind solving systems of linear equations of the form Ax=b where
  • A is an m×n matrix of coefficients for m equations in n unknowns
  • x is an n×1 vector unknowns, x1,x2,xn
  • b is an m×1 vector of constants, the “right-hand sides” of the equations
The general conditions for solutions are:
  • the equations are consistent (solutions exist) if r(A|b)=r(A)
    • the solution is unique if r(A|b)=r(A)=n
    • the solution is underdetermined if r(A|b)=r(A)<n
  • the equations are inconsistent (no solutions) if r(A|b)>r(A)
We use c( R(A), R(cbind(A,b)) ) to show the ranks, and all.equal( R(A), R(cbind(A,b)) ) to test for consistency.
library(matlib)   # use the package

Equations in two unknowns

Each equation in two unknowns corresponds to a line in 2D space. The equations have a unique solution if all lines intersect in a point.

Two consistent equations

A <- matrix(c(1, 2, -1, 2), 2, 2)
b <- c(2,1)
showEqn(A, b)
## 1*x1 - 1*x2  =  2 
## 2*x1 + 2*x2  =  1
c( R(A), R(cbind(A,b)) )          # show ranks
## [1] 2 2
all.equal( R(A), R(cbind(A,b)) )  # consistent?
## [1] TRUE
Plot them:
plotEqn(A,b)
##   x1   - x2  =  2 
## 2*x1 + 2*x2  =  1

Three consistent equations

For three (or more) equations in two unknowns, r(A)2, because r(A)min(m,n). The equations will be consistent if r(A)=r(A|b), which means that whatever linear relations exist among the rows of A are the same as those among the elements of b.
A <- matrix(c(1,2,3, -1, 2, 1), 3, 2)
b <- c(2,1,3)
showEqn(A, b)
## 1*x1 - 1*x2  =  2 
## 2*x1 + 2*x2  =  1 
## 3*x1 + 1*x2  =  3
c( R(A), R(cbind(A,b)) )          # show ranks
## [1] 2 2
all.equal( R(A), R(cbind(A,b)) )  # consistent?
## [1] TRUE
Plot them:
plotEqn(A,b)
##   x1   - x2  =  2 
## 2*x1 + 2*x2  =  1 
## 3*x1   + x2  =  3

Three inconsistent equations

A <- matrix(c(1,2,3, -1, 2, 1), 3, 2)
b <- c(2,1,6)
showEqn(A, b)
## 1*x1 - 1*x2  =  2 
## 2*x1 + 2*x2  =  1 
## 3*x1 + 1*x2  =  6
c( R(A), R(cbind(A,b)) )          # show ranks
## [1] 2 3
all.equal( R(A), R(cbind(A,b)) )  # consistent?
## [1] "Mean relative difference: 0.5"
An approximate solution is sometimes available using a generalized inverse.
x <- MASS::ginv(A) %*% b
x
##      [,1]
## [1,]    2
## [2,]   -1
Plot the equations. You can see that each pair of equations has a solution, but all three do not have a common, consistent solution.
plotEqn(A,b, xlim=c(-2, 4))
##   x1   - x2  =  2 
## 2*x1 + 2*x2  =  1 
## 3*x1   + x2  =  6
points(x[1], x[2], pch=15)

Equations in three unknowns

Each equation in three unknowns corresponds to a plane in 3D space. The equations have a unique solution if all planes intersect in a point.

Three consistent equations

A <- matrix(c(2, 1, -1,
             -3, -1, 2,
             -2,  1, 2), 3, 3, byrow=TRUE)
colnames(A) <- paste0('x', 1:3)
b <- c(8, -11, -3)
showEqn(A, b)
##  2*x1 + 1*x2 - 1*x3  =    8 
## -3*x1 - 1*x2 + 2*x3  =  -11 
## -2*x1 + 1*x2 + 2*x3  =   -3
Are the equations consistent?
c( R(A), R(cbind(A,b)) )          # show ranks
## [1] 3 3
all.equal( R(A), R(cbind(A,b)) )  # consistent?
## [1] TRUE
Solve for x.
solve(A, b)
## x1 x2 x3 
##  2  3 -1
solve(A) %*% b
##    [,1]
## x1    2
## x2    3
## x3   -1
inv(A) %*% b
##      [,1]
## [1,]    2
## [2,]    3
## [3,]   -1
Another way to see the solution is to reduce A|b to echelon form. The result is I|A1b, with the solution in the last column.
echelon(A, b)
##      x1 x2 x3   
## [1,]  1  0  0  2
## [2,]  0  1  0  3
## [3,]  0  0  1 -1
echelon(A, b, verbose=TRUE, fractions=TRUE)
## 
## Initial matrix:
##      x1  x2  x3     
## [1,]   2   1  -1   8
## [2,]  -3  -1   2 -11
## [3,]  -2   1   2  -3
## 
## row: 1 
## 
##  exchange rows 1 and 2 
##      x1  x2  x3     
## [1,]  -3  -1   2 -11
## [2,]   2   1  -1   8
## [3,]  -2   1   2  -3
## 
##  multiply row 1 by -1/3 
##      x1   x2   x3       
## [1,]    1  1/3 -2/3 11/3
## [2,]    2    1   -1    8
## [3,]   -2    1    2   -3
## 
##  multiply row 1 by 2 and subtract from row 2 
##      x1   x2   x3       
## [1,]    1  1/3 -2/3 11/3
## [2,]    0  1/3  1/3  2/3
## [3,]   -2    1    2   -3
## 
##  multiply row 1 by 2 and add to row 3 
##      x1   x2   x3       
## [1,]    1  1/3 -2/3 11/3
## [2,]    0  1/3  1/3  2/3
## [3,]    0  5/3  2/3 13/3
## 
## row: 2 
## 
##  exchange rows 2 and 3 
##      x1   x2   x3       
## [1,]    1  1/3 -2/3 11/3
## [2,]    0  5/3  2/3 13/3
## [3,]    0  1/3  1/3  2/3
## 
##  multiply row 2 by 3/5 
##      x1   x2   x3       
## [1,]    1  1/3 -2/3 11/3
## [2,]    0    1  2/5 13/5
## [3,]    0  1/3  1/3  2/3
## 
##  multiply row 2 by 1/3 and subtract from row 1 
##      x1   x2   x3       
## [1,]    1    0 -4/5 14/5
## [2,]    0    1  2/5 13/5
## [3,]    0  1/3  1/3  2/3
## 
##  multiply row 2 by 1/3 and subtract from row 3 
##      x1   x2   x3       
## [1,]    1    0 -4/5 14/5
## [2,]    0    1  2/5 13/5
## [3,]    0    0  1/5 -1/5
## 
## row: 3 
## 
##  multiply row 3 by 5 
##      x1   x2   x3       
## [1,]    1    0 -4/5 14/5
## [2,]    0    1  2/5 13/5
## [3,]    0    0    1   -1
## 
##  multiply row 3 by 4/5 and add to row 1 
##      x1   x2   x3       
## [1,]    1    0    0    2
## [2,]    0    1  2/5 13/5
## [3,]    0    0    1   -1
## 
##  multiply row 3 by 2/5 and subtract from row 2 
##      x1 x2 x3   
## [1,]  1  0  0  2
## [2,]  0  1  0  3
## [3,]  0  0  1 -1
Plot them. plotEqn3d uses rgl for 3D graphics. If you rotate the figure, you’ll see an orientation where all three planes intersect at the solution point, x=(2,3,1)
plotEqn3d(A,b, xlim=c(0,4), ylim=c(0,4))

Three inconsistent equations

A <- matrix(c(1,  3, 1,
              1, -2, -2,
              2,  1, -1), 3, 3, byrow=TRUE)
colnames(A) <- paste0('x', 1:3)
b <- c(2, 3, 6)
showEqn(A, b)
## 1*x1 + 3*x2 + 1*x3  =  2 
## 1*x1 - 2*x2 - 2*x3  =  3 
## 2*x1 + 1*x2 - 1*x3  =  6
Are the equations consistent? No.
c( R(A), R(cbind(A,b)) )          # show ranks
## [1] 2 3
all.equal( R(A), R(cbind(A,b)) )  # consistent?
## [1] "Mean relative difference: 0.5"

2016年11月6日星期日

Placing Text in Plots

get( getOption( "device" ) )()

##  Create the plot. This sets the "usr" coordinates.

plot( 1:10 )

##  par( "usr" ) returns a vector containing xleft, xright, ybottom, ytop.

usr <- par( "usr" )

##  Place the text. Note the use of the adj parameter.
##  Left top corner:

text( usr[ 1 ], usr[ 4 ], "left top",     adj = c( 0, 1 ), col = "blue" )

##  Left bottom corner:

text( usr[ 1 ], usr[ 3 ], "left bottom",  adj = c( 0, 0 ), col = "blue" )

##  Right top corner:

text( usr[ 2 ], usr[ 4 ], "right top",    adj = c( 1, 1 ), col = "blue" )

##  Right bottom corner:

text( usr[ 2 ], usr[ 3 ], "right bottom", adj = c( 1, 0 ), col = "blue" )
http://sphaerula.com/legacy/R/placingTextInPlots.html