knitr::opts_chunk$set( warning = FALSE, message = FALSE ) options(digits=4) par(mar=c(5,4,1,1)+.1)
This example shows two classical ways to find the determinant, $\det(A)$ of a square matrix. They each work by reducing the problem to a series of smaller ones which can be more easily calculated.
library(matlib)
det()
by cofactor expansionSet up a $3 \times 3$ matrix, and find its determinant (so we know what the answer should be).
A <- matrix(c(4, 2, 1, 5, 6, 7, 1, 0, 3), nrow=3, byrow=TRUE) det(A)
The cofactor $A_{i,j}$ of element $a_{i,j}$ is the signed determinant of what is left when row i, column j of the matrix $A$ are deleted. NB: In R, negative subscripts delete rows or columns.
cat(cofactor(A, 1, 1), " == ", 1 * det( (A[-1, -1]), "\n" )) cat(cofactor(A, 1, 2), " == ", -1 * det( (A[-1, -2]), "\n" )) cat(cofactor(A, 1, 3), " == ", 1 * det( (A[-1, -3]), "\n" ))
In symbols: $\det(A) = a_{1,1} * A_{1,1} + a_{1,2} * A_{1,2} + a_{1,3} * A_{1,3}$
rowCofactors()
is a convenience function, that calculates these all together
rowCofactors(A, 1)
Voila: Multiply row 1 times the cofactors of its elements. NB: In R, this multiplication gives a $1 \times 1$ matrix.
A[1,] %*% rowCofactors(A, 1) all.equal( det(A), c(A[1,] %*% rowCofactors(A, 1)) )
det()
by Gaussian elimination (pivoting)This example follows Green and Carroll, Table 2.2.
Start with a 4 x 4 matrix, $M$, and save det(M)
.
M <- matrix(c(2, 3, 1, 2, 4, 2, 3, 4, 1, 4, 2, 2, 3, 1, 0, 1), nrow=4, ncol=4, byrow=TRUE) (dsave <- det(M)) # ### 'pivot' on the leading diagonal element, M[1,1]:
det()
will be the product of the 'pivots', the leading diagonal elements.
This step reduces row 1 and column 1 to 0, so it may be discarded.
NB: In R, dropping a row/column can change a matrix to a vector, so we use
drop = FALSE
inside the subscript.
(d <- M[1,1]) #-- Reduce row 1, col 1 to 0 (M[1,] <- M[1,, drop=FALSE] / M[1, 1]) (M <- M - M[,1] %*% M[1,, drop=FALSE]) #-- Drop first row and column M <- M[-1, -1] #-- Accumulate the product of pivots d <- d * M[1, 1]
(M[1,] <- M[1,, drop=FALSE] / M[1,1]) (M <- M - M[,1] %*% M[1,, drop=FALSE]) M <- M[-1, -1] d = d * M[1, 1]
(M[1,] <- M[1,, drop=FALSE] / M[1,1]) (M <- M - M[,1] %*% M[1,, drop=FALSE]) M <- M[-1, -1, drop=FALSE] d <- d * M[1, 1] # did we get it right? all.equal(d, dsave)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.