2. Evaluation of determinants

Michael Friendly

2024年11月24日

Source: vignettes/a2-det-ex2.Rmd

a2-det-ex2.Rmd

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 )

1. Calculate `det()` by cofactor expansion

Set 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)

## [1] 50

Find cofactors of row 1 elements

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" ))

## 18 == 18

 cat (cofactor (A, 1, 2), " == ", -1 * det ( (A[-1, -2]), "\n" ))

## -8 == -8

 cat (cofactor (A, 1, 3), " == ", 1 * det ( (A[-1, -3]), "\n" ))

## -6 == -6

det() = product of row with cofactors

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)

## [1] 18 -8 -6

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)

## [,1]
## [1,] 50

 all.equal ( det (A), c (A[1,] %*%  rowCofactors (A, 1)) )

## [1] TRUE

2. Finding `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))

## [1] 15

# ### '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])

## [1] 2

 #-- Reduce row 1, col 1 to 0
 (M[1,] <- M[1,, drop=FALSE] / M[1, 1])

## [,1] [,2] [,3] [,4]
## [1,] 1 1.5 0.5 1

 (M <- M - M[,1] %*%  M[1,, drop=FALSE])

## [,1] [,2] [,3] [,4]
## [1,] 0 0.0 0.0 0
## [2,] 0 -4.0 1.0 0
## [3,] 0 2.5 1.5 1
## [4,] 0 -3.5 -1.5 -2

 #-- Drop first row and column
 M <- M[-1, -1]

 #-- Accumulate the product of pivots
 d <- d * M[1, 1]

Repeat, reducing new row, col 1 to 0

 (M[1,] <- M[1,, drop=FALSE] / M[1,1])

## [,1] [,2] [,3]
## [1,] 1 -0.25 0

 (M <- M - M[,1] %*%  M[1,, drop=FALSE])

## [,1] [,2] [,3]
## [1,] 0 0.000 0
## [2,] 0 2.125 1
## [3,] 0 -2.375 -2

 M <- M[-1, -1]
 d = d * M[1, 1]

Repeat once more. d = det(M)

 (M[1,] <- M[1,, drop=FALSE] / M[1,1])

## [,1] [,2]
## [1,] 1 0.4706

 (M <- M - M[,1] %*%  M[1,, drop=FALSE])

## [,1] [,2]
## [1,] 0 0.0000
## [2,] 0 -0.8824

 M <- M[-1, -1, drop=FALSE]
 d <- d * M[1, 1]

 # did we get it right?
 all.equal (d, dsave)

## [1] TRUE