Evaluation of determinants"

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)

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)

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" ))
  cat(cofactor(A, 1, 2),  "  ==  ", -1 * det( (A[-1, -2]), "\n" ))
  cat(cofactor(A, 1, 3),  "  ==  ",  1 * det( (A[-1, -3]), "\n" ))

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)

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

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

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

Repeat, reducing new row, col 1 to 0

  (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]

Repeat once more. d = det(M)

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


Try the matlib package in your browser

Any scripts or data that you put into this service are public.

matlib documentation built on April 4, 2018, 5:03 p.m.