ASSIGNMENT 2 IS 605 FUNDAMENTALS OF COMPUTATIONAL MATHEMATICS - 2015 1. Problem Set 1 (1) Show that \(A^TA\neq AA^T\) in general. (Proof and demonstration.) (2) For a special type of square matrix A, we get \(A^TA=AA^T\). Under what conditions could this be true? (Hint: The Identity matrix I is an example of such a matrix).
Let \(A=\left(a_{ij}\right)\in\mathbb{R}^{m\times n}\) and \(A^T=\left(a_{ji}\right)\in\mathbb{R}^{n\times m}\). To compute the product \(A^TA\), the entry in the product at row \(i\) and column \(j\) is given by the dot product of the \(i\)-th row of \(A^T\) and the \(j\)-th column of \(A\) = \(\sum_{k<m} a_{ki}a_{kj}\), so \[A^TA = \left(\sum_{k<m} a_{ki}a_{kj}\right)\in\mathbb{R}^{n\times n}.\] Similarly, the dot product of the \(i\)-th row of \(A\) and the \(j\)-th column of \(A^T\) is \(\sum_{k<n}a_{ik}a_{jk}\) so \[AA^T=\left(\sum_{k<n}a_{ik}a_{jk}\right)\in\mathbb{R}^{m\times m}.\] It is clear, then, that for \(n\neq m\), the dimensions of \(A^TA\) and \(AA^T\) do not even match, so they cannot be equal. When \(n=m\), \(A^TA=AA^T\) if and only if for all \(i,j<n\), \[\begin{equation}\label{condition} \sum_{k<n}a_{ki}a_{kj}=\sum_{k<n}a_{ik}a_{jk}\tag{*} \end{equation}\]
As a concrete example, consider the matrix, \[A= \begin{pmatrix} 1 & 1\\ 0 & 0 \end{pmatrix}\in\mathbb{R}^{2\times 2}.\] It is easy to see that for \(i=0\) and \(j=0\), \[\begin{align*} \sum_{k<2}a_{k0}a_{k0} &= a_{00}a_{00}+a_{10}a_{10} \\ &= 1\cdot 1 + 0\cdot 0 \\ &= 1 \end{align*}\] however, \[\begin{align*} \sum_{k<2}a_{0k}a_{0k} &= a_{00}a_{00}+a_{01}a_{01} \\ &= 1\cdot 1 + 1\cdot 1 \\ &= 2 \end{align*}\] And we see that \[A^TA= \begin{pmatrix} 1 & 1\\ 1 & 1 \end{pmatrix}\] \[AA^T= \begin{pmatrix} 2 & 0\\ 0 & 0 \end{pmatrix}\]
For a special type of square matrix \(A\), we get \(A^TA=AA^T\). Under what conditions could this be true? (Hint: The Identity matrix \(I\) is an example of such a matrix.)\
As seen in problem (1), for an \(n\times n\) matrix \(A=\left(a_{ij}\right)\), \(A^TA=AA^T\) if and only if for all \(i,j<n\), \(A\) satisfies equation .
One particular class of such matricies are symmetric matricies, that is, all square matricies \(A\) such that \(A=A^T\). For such matricies, it is easy to see that \[A^TA=AA=AA^T.\] Given a square matrix \(A=\left(a_{ij}\right)\in\mathbb{R}^{n\times n}\), since \(A^T=\left(a_{ji}\right)\), symmetry is equivalent to the statement, for all \(i,j<n\), \(a_{ij}=a_{ji}\). From that it is easy to see that for such a matrix, \[\sum_{k<n}a_{ki}a_{kj}=\sum_{k<n}a_{ik}a_{jk}\] and thus satisfies equation .
However, there are matricies that satisfy equation but are not symmetric, e.g. \[A= \begin{pmatrix} 1 & 1 & 0\\ 0 & 1 & 1\\ 1 & 0 & 1 \end{pmatrix}\]
It is easy to see that \(A\) is symmetric (although its rows can be shuffled to make it symmetric), however, \[A^TA= \begin{pmatrix} 2 & 1 & 1\\ 1 & 2 & 1\\ 1 & 1 & 2 \end{pmatrix}=AA^T.\]
My solution to this problem borrows very heavily from the matrixcalc package’s function ‘lu.decomposition’ lu.decomposition
This solution uses Doolittle decomposition method.
This method LU decomposition is assumed to exist. We start with the form of L and U. We then solve for the elements in L and U from the equations that result from the multiplications necessary for A = LU
lu.decomp <- function(x){
#count number of rows in matrix
n <- nrow( x )
#create matrix of zeroes as starting point for both the lower and upper triangular matrices
L <- matrix( 0, nrow=n, ncol=n )
U <- matrix( 0, nrow=n, ncol=n )
#create identity matrix
diag( L ) <- rep( 1, n )
#iterate through each row
for ( i in 1:n ) {
ip1 <- i + 1
im1 <- i - 1
#iterate through each column
for ( j in 1:n ) {
#set the element in U to the value in the original matrix
U[i,j] <- x[i,j]
#summation of L(i,j) * U(i,j) and evaluation of U(i,j)
if ( im1 > 0 ) {
for ( k in 1:im1 ) {
U[i,j] <- U[i,j] - L[i,k] * U[k,j]
}
}
}
if ( ip1 <= n ) {
for ( j in ip1:n ) {
L[j,i] <- x[j,i]
if ( im1 > 0 ) {
#summation of L(i,j) * U(i,j)
for ( k in 1:im1 ) {
L[j,i] <- L[j,i] - L[j,k] * U[k,i]
}
}
if ( U[i,i] == 0 )
stop( "argument x is a singular matrix" )
##Evaluation of L(i,j)
L[j,i] <- L[j,i] / U[i,i]
}
}
}
result <- list( L=L, U=U )
return( result )
}
#Create a test matrix
mat1.data <- c(1,2,3,4,5,6,7,8,9)
mat1 <- matrix(mat1.data,nrow=3,ncol=3,byrow=TRUE)
mat1## [,1] [,2] [,3]
## [1,] 1 2 3
## [2,] 4 5 6
## [3,] 7 8 9lu.decomp(mat1) ## $L
## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 4 1 0
## [3,] 7 2 1
##
## $U
## [,1] [,2] [,3]
## [1,] 1 2 3
## [2,] 0 -3 -6
## [3,] 0 0 0