Show proof that \(A^TA\) \(\neq\) \(AA^T\) in general.
A = \(\begin{bmatrix} a & b &c \\ d & e & f \end{bmatrix}\)
\(A^T\) = \(\begin{bmatrix} a & d \\ b & e \\ c & f \end{bmatrix}\)
First, I will show \(A^TA\).
\(\begin{bmatrix} a & d \\ b & e \\ c & f \end{bmatrix}\) * \(\begin{bmatrix} a & b &c \\ d & e & f \end{bmatrix}\)
\(\begin{bmatrix} aa + dd & ab + de & ac + df \\ ba + ed & bb + ee & bc + ef \\ ca + fd & cb + fe & cc + ff \end{bmatrix}\)
Next, I will show \(AA^T\).
\(\begin{bmatrix} a & b &c \\ d & e & f \end{bmatrix}\) * \(\begin{bmatrix} a & d \\ b & e \\ c & f \end{bmatrix}\)
\(\begin{bmatrix} aa + bb + cc & ad + be + cf \\ da + eb + fc & dd + ee + ff \end{bmatrix}\)
Now, I will show them side by side.
\(A^TA\) \(AA^T\)
\(\begin{bmatrix} aa + dd & ab + de & ac + df \\ ba + ed & bb + ee & bc + ef \\ ca + fd & cb + fe & cc + ff \end{bmatrix}\) . . . \(\begin{bmatrix} aa + bb + cc & ad + be + cf \\ da + eb + fc & dd + ee + ff \end{bmatrix}\)
Therefore, \(A^TA\) will not always equal \(AA^T\).
An example of this is shown below.
A = \(\begin{bmatrix} 2 & 4 & 8 \\ 3 & 1 & 7 \end{bmatrix}\)
\(A^T\) = \(\begin{bmatrix} 2 & 3 \\ 4 & 1 \\ 8 & 7 \end{bmatrix}\)
First, I will show \(A^TA\).
\(\begin{bmatrix} 2 & 3 \\ 4 & 1 \\ 8 & 7 \end{bmatrix}\) * \(\begin{bmatrix} 2 & 4 & 8 \\ 3 & 1 & 7 \end{bmatrix}\)
\(\begin{bmatrix} 2*2 + 3*3 & 2*4 + 3*1 & 2*8 + 3*7 \\ 4*2 + 1*3 & 4*4 + 1*1 & 4*8 + 1*7 \\ 8*2 + 7*3 & 8*4 + 7*1 & 8*8 + 8*7 \end{bmatrix}\)
\(\begin{bmatrix} 13 & 11 & 37 \\ 11 & 17 & 39 \\ 37 & 39 & 120 \end{bmatrix}\)
Next, I will show \(AA^T\).
\(\begin{bmatrix} 2 & 4 & 8 \\ 3 & 1 & 7 \end{bmatrix}\) * \(\begin{bmatrix} 2 & 3 \\ 4 & 1 \\ 8 & 7 \end{bmatrix}\)
\(\begin{bmatrix} 2*2 + 4*4 + 8*8 & 2*3 + 4*1 + 8*7 \\ 3*2 + 1*4 + 7*8 & 3*3 + 1*1 + 7*7 \end{bmatrix}\)
\(\begin{bmatrix} 84 & 66 \\ 66 & 59 \end{bmatrix}\)
Now, I will show them side by side.
\(A^TA\) \(AA^T\)
\(\begin{bmatrix} 13 & 11 & 37 \\ 11 & 17 & 39 \\ 37 & 39 & 120 \end{bmatrix}\) . . . \(\begin{bmatrix} 84 & 66 \\ 66 & 59 \end{bmatrix}\)
You can see with this example that \(A^TA\) will not always equal \(AA^T\).
For a special type of square matrix A, we get \(A^TA\) = \(AA^T\) . Under what conditions could this be true?
The following bullet points show two requirements for this statement to be true:
A matrix that is symmetrical will make this statement true. The reason for this is that \(A^T\) will be equivalent to \(A\) when the matrix is symmetrical. For example, a matrix of \(\begin{bmatrix} a & b & c \\ b & d & e \\ c & e & f \end{bmatrix}\) will have a transpose matrix of \(\begin{bmatrix} a & b & c \\ b & d & e \\ c & e & f \end{bmatrix}\). When \(A^T\) is the same as \(A\), then \(A^TA\) = \(AA^T\) becomes \(AA = AA\), which is true.
The matrix will need to be a square matrix. If the matrix is not a square matrix, then each side of the equation will have different numbers of rows and columns. For example a matrix of 3 x 2 will have a transpose of a matrix of 2 x 3. When they are multiplied together in order, the resulting matrix will be 3 x 3. When they are multiplied the opposite order, the resulting matrix will be 2 x 2. Different numbers of rows and columns means different matrices. Therefore, the matrix will need to be a square matrix.
Write an R function to factorize a square matrix A into LU or LDU, whichever you prefer.
In this example, I created a function to factorize a square matrix into LU. After creating the function, I show its ability to work with many examples.
My function in the first R chunk only works for matrices that have the dimensions 3 x 3 or 4 x 4. I did not add a section that would work for matrices with a dimension of 5 x 5 because the first two types of square matrices contained plenty of code to demonstrate my understanding.
library(matrixcalc)
lufactorization <- function(myMatrix)
{
l1 <- 0
l2 <- 0
l3 <- 0
p <- matrix(c(1,0,0,0,1,0,0,0,1))
useP <- FALSE
saveForLater <- myMatrix
skip2 <- FALSE
skipRest <- FALSE
if (is.square.matrix(myMatrix) == TRUE)
{
if (nrow(myMatrix) == 3)
{
d <- myMatrix[2,1]
if ( ((myMatrix[1,1] > 0)&(d > 0)) || ((myMatrix[1,1] < 0)&(d < 0)) )
{
myMatrix <- matrix(c(myMatrix[1,1], myMatrix[1,2], myMatrix[1,3], myMatrix[2,1] - (myMatrix[1,1] * d / myMatrix[1,1]), myMatrix[2,2] - (myMatrix[1,2] * d / myMatrix[1,1]), myMatrix[2,3] - (myMatrix[1,3] * d / myMatrix[1,1]), myMatrix[3,1], myMatrix[3,2], myMatrix[3,3]), nrow = 3, ncol = 3, byrow = TRUE)
l1 <- d / myMatrix[1,1]
}
else if ((myMatrix[1,1] != 0) & (d != 0))
{
myMatrix <- matrix(c(myMatrix[1,1], myMatrix[1,2], myMatrix[1,3], myMatrix[2,1] - (myMatrix[1,1] * d / myMatrix[1,1]), myMatrix[2,2] - (myMatrix[1,2] * d / myMatrix[1,1]), myMatrix[2,3] - (myMatrix[1,3] * d / myMatrix[1,1]), myMatrix[3,1], myMatrix[3,2], myMatrix[3,3]), nrow = 3, ncol = 3, byrow = TRUE)
l1 <- d / myMatrix[1,1]
}
else
{
if (myMatrix[2,1] != 0)
{
myMatrix <- matrix(c(myMatrix[2,1], myMatrix[2,2], myMatrix[2,3], myMatrix[1,1], myMatrix[1,2], myMatrix[1,3], myMatrix[3,1], myMatrix[3,2], myMatrix[3,3]), nrow = 3, ncol = 3, byrow = TRUE)
l1 <- 0
useP <- TRUE
p <- matrix(c(0,1,0,1,0,0,0,0,1), nrow = 3, ncol = 3, byrow = TRUE)
}
else if (myMatrix[3,1] != 0)
{
myMatrix <- matrix(c(myMatrix[3,1], myMatrix[3,2], myMatrix[3,3], myMatrix[2,1], myMatrix[2,2], myMatrix[2,3], myMatrix[1,1], myMatrix[1,2], myMatrix[1,3]), nrow = 3, ncol = 3, byrow = TRUE)
l1 <- 0
l2 <- 0
useP <- TRUE
p <- matrix(c(0,0,1,0,1,0,1,0,0))
skip2 <- TRUE
}
else
{
print(" ")
print("Your original matrix")
print(saveForLater)
threeIMatrix <- matrix(c(1, 0, 0, 0, 1, 0, 0, 0, 1), nrow = 3, ncol = 3, byrow = TRUE)
print("Your matrix L is")
print(threeIMatrix)
print("Your matrix U is")
print(myMatrix)
print("This cannot be decomposed any further.")
return()
}
}
if (skip2 != TRUE)
{
g <- myMatrix[3,1]
if ( ((myMatrix[1,1] > 0)&(g > 0)) || ((myMatrix[1,1] < 0)&(g < 0)) )
{
myMatrix <- matrix(c(myMatrix[1,1], myMatrix[1,2], myMatrix[1,3], myMatrix[2,1], myMatrix[2,2], myMatrix[2,3], myMatrix[3,1] - (myMatrix[1,1] * g / myMatrix[1,1]), myMatrix[3,2] - (myMatrix[1,2] * g / myMatrix[1,1]), myMatrix[3,3] - (myMatrix[1,3] * g / myMatrix[1,1])), nrow = 3, ncol = 3, byrow = TRUE)
l2 <- g / myMatrix[1,1]
}
else if ((myMatrix[1,1] != 0) & (g != 0))
{
myMatrix <- matrix(c(myMatrix[1,1], myMatrix[1,2], myMatrix[1,3], myMatrix[2,1], myMatrix[2,2], myMatrix[2,3], myMatrix[3,1] - (myMatrix[1,1] * g / myMatrix[1,1]), myMatrix[3,2] - (myMatrix[1,2] * g / myMatrix[1,1]), myMatrix[3,3] - (myMatrix[1,3] * g / myMatrix[1,1])), nrow = 3, ncol = 3, byrow = TRUE)
l2 <- g / myMatrix[1,1]
}
else
{
l2 <- 0
}
}
h <- myMatrix[3,2]
if ( ((myMatrix[2,2] > 0)&(h > 0)) || ((myMatrix[2,2] < 0)&(h < 0)) )
{
myMatrix <- matrix(c(myMatrix[1,1], myMatrix[1,2], myMatrix[1,3], myMatrix[2,1], myMatrix[2,2], myMatrix[2,3], myMatrix[3,1], myMatrix[3,2] - (myMatrix[2,2] * h / myMatrix[2,2]), myMatrix[3,3] - (myMatrix[2,3] * h / myMatrix[2,2])), nrow = 3, ncol = 3, byrow = TRUE)
l3 <- h / myMatrix[2,2]
}
else if ((myMatrix[2,2] != 0) & (h != 0))
{
myMatrix <- matrix(c(myMatrix[1,1], myMatrix[1,2], myMatrix[1,3], myMatrix[2,1], myMatrix[2,2], myMatrix[2,3], myMatrix[3,1], myMatrix[3,2] - (myMatrix[2,2] * h / myMatrix[2,2]), myMatrix[3,3] - (myMatrix[2,3] * h / myMatrix[2,2])), nrow = 3, ncol = 3, byrow = TRUE)
l3 <- h / myMatrix[2,2] * -1
}
else
{
l3 <- 0
}
print("Your original matrix")
print(saveForLater)
print("Your matrix L is")
lMatrix <- matrix(c(1,0,0,l1,1,0,l2,l3,0), nrow = 3, ncol = 3, byrow = TRUE)
print(lMatrix)
print("Your matrix U is")
print(myMatrix)
if (useP == TRUE)
{
print("Your matrix P is")
print(p)
}
}
else if (nrow(myMatrix) == 4)
{
e <- myMatrix[2,1]
if ( ((myMatrix[1,1] > 0)&(e > 0)) || ((myMatrix[1,1] < 0)&(e < 0)) )
{
myMatrix <- matrix(c(myMatrix[1,1], myMatrix[1,2], myMatrix[1,3], myMatrix[1,4], myMatrix[2,1] - (myMatrix[1,1] * e / myMatrix[1,1]), myMatrix[2,2] - (myMatrix[1,2] * e / myMatrix[1,1]), myMatrix[2,3] - (myMatrix[1,3] * e / myMatrix[1,1]), myMatrix[2,4] -(myMatrix[1,4] * e / myMatrix[1,1]), myMatrix[3,1], myMatrix[3,2], myMatrix[3,3], myMatrix[3,4], myMatrix[4,1], myMatrix[4,2], myMatrix[4,3], myMatrix[4,4]), nrow = 4, ncol = 4, byrow = TRUE)
l1 <- e / myMatrix[1,1]
}
else if ((myMatrix[1,1] != 0) & (e != 0))
{
myMatrix <- matrix(c(myMatrix[1,1], myMatrix[1,2], myMatrix[1,3], myMatrix[1,4], myMatrix[2,1] - (myMatrix[1,1] * e / myMatrix[1,1]), myMatrix[2,2] - (myMatrix[1,2] * e / myMatrix[1,1]), myMatrix[2,3] - (myMatrix[1,3] * e / myMatrix[1,1]), myMatrix[2,4] - (myMatrix[1,4] * e / myMatrix[1,1]), myMatrix[3,1], myMatrix[3,2], myMatrix[3,3], myMatrix[3,4], myMatrix[4,1], myMatrix[4,2], myMatrix[4,3], myMatrix[4,4]), nrow = 4, ncol = 4, byrow = TRUE)
l1 <- e / myMatrix[1,1]
}
else
{
if (myMatrix[2,1] != 0)
{
myMatrix <- matrix(c(myMatrix[2,1], myMatrix[2,2], myMatrix[2,3], myMatrix[2,4], myMatrix[1,1], myMatrix[1,2], myMatrix[1,3], myMatrix[1,4], myMatrix[3,1], myMatrix[3,2], myMatrix[3,3], myMatrix[3,4], myMatrix[4,1], myMatrix[4,2], myMatrix[4,3], myMatrix[4,4]), nrow = 4, ncol = 4, byrow = TRUE)
l1 <- 0
useP <- TRUE
p <- matrix(c(0,1,0,0,1,0,0,0,0,0,1,0,0,0,0,1), nrow = 4, ncol = 4, byrow = TRUE)
}
else if (myMatrix[3,1] != 0)
{
myMatrix <- matrix(c(myMatrix[3,1], myMatrix[3,2], myMatrix[3,3], myMatrix[3,4], myMatrix[2,1], myMatrix[2,2], myMatrix[2,3], myMatrix[2,4], myMatrix[3,1], myMatrix[1,2], myMatrix[1,3], myMatrix[1,4], myMatrix[1,1], myMatrix[4,2], myMatrix[4,3], myMatrix[4,4]), nrow = 4, ncol = 4, byrow = TRUE)
l1 <- 0
useP <- TRUE
p <- matrix(c(0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,1), nrow = 4, ncol = 4, byrow = TRUE)
}
else
{
print("Your original matrix")
print(saveForLater)
fourIMatrix <- matrix(c(1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1), nrow = 4, ncol = 4, byrow = TRUE)
print("Your matrix L is")
print(fourIMatrix)
print("Your matrix U is")
print(myMatrix)
print("This cannot be decomposed any further.")
return()
}
}
i <- myMatrix[3,1]
if ( ((myMatrix[1,1] > 0)&(i > 0)) || ((myMatrix[1,1] < 0)&(i < 0)) )
{
myMatrix <- matrix(c(myMatrix[1,1], myMatrix[1,2], myMatrix[1,3], myMatrix[1,4], myMatrix[2,1], myMatrix[2,2], myMatrix[2,3], myMatrix[2,4], myMatrix[3,1] - (myMatrix[1,1] * i / myMatrix[1,1]), myMatrix[3,2] - (myMatrix[1,2] * i / myMatrix[1,1]), myMatrix[3,3] - (myMatrix[1,3] * i / myMatrix[1,1]), myMatrix[3,4] - (myMatrix[1,4] * i / myMatrix[1,1]), myMatrix[4,]), nrow = 4, ncol = 4, byrow = TRUE)
l2 <- i / myMatrix[1,1]
}
else if ((myMatrix[1,1] != 0) & (i != 0))
{
myMatrix <- matrix(c(myMatrix[1,1], myMatrix[1,2], myMatrix[1,3], myMatrix[1,4], myMatrix[2,1], myMatrix[2,2], myMatrix[2,3], myMatrix[2,4], myMatrix[3,1] - (myMatrix[1,1] * i / myMatrix[1,1]), myMatrix[3,2] - (myMatrix[1,2] * i / myMatrix[1,1]), myMatrix[3,3] - (myMatrix[1,3] * i / myMatrix[1,1]), myMatrix[3,4] - (myMatrix[1,4] * i / myMatrix[1,1]), myMatrix[4,]), nrow = 4, ncol = 4, byrow = TRUE)
l2 <- i / myMatrix[1,1]
}
else
{
l2 <- 0
}
j <- myMatrix[3,2]
if ( ((myMatrix[2,2] > 0)&(j > 0)) || ((myMatrix[2,2] < 0)&(j < 0)) )
{
myMatrix <- matrix(c(myMatrix[1,1], myMatrix[1,2], myMatrix[1,3], myMatrix[1,4], myMatrix[2,1], myMatrix[2,2], myMatrix[2,3], myMatrix[2,4], myMatrix[3,1], myMatrix[3,2] - (myMatrix[2,2] * j / myMatrix[2,2]), myMatrix[3,3] - (myMatrix[2,3] * j / myMatrix[2,2]), myMatrix[3,4] - myMatrix[2,4] * j / myMatrix[2,2], myMatrix[4,]), nrow = 4, ncol = 4, byrow = TRUE)
l3 <- j / myMatrix[2,2]
}
else if ((myMatrix[2,2] != 0) & (j != 0))
{
myMatrix <- matrix(c(myMatrix[1,1], myMatrix[1,2], myMatrix[1,3], myMatrix[1,4], myMatrix[2,1], myMatrix[2,2], myMatrix[2,3], myMatrix[2,4], myMatrix[3,1], myMatrix[3,2] - (myMatrix[2,2] * j / myMatrix[2,2]), myMatrix[3,3] - (myMatrix[2,3] * j / myMatrix[2,2]), myMatrix[3,4] - myMatrix[2,4] * j / myMatrix[2,2], myMatrix[4,]), nrow = 4, ncol = 4, byrow = TRUE)
l3 <- j / myMatrix[2,2]
}
else
{
l3 <- 0
}
m <- myMatrix[4,1]
if ( ((myMatrix[1,1] > 0)&(m > 0)) || ((myMatrix[1,1] < 0)&(m < 0)) )
{
myMatrix <- matrix(c(myMatrix[1,1], myMatrix[1,2], myMatrix[1,3], myMatrix[1,4], myMatrix[2,1], myMatrix[2,2], myMatrix[2,3], myMatrix[2,4], myMatrix[3,], myMatrix[4,1] - (myMatrix[1,1] * m / myMatrix[1,1]), myMatrix[4,2] - (myMatrix[1,2] * m / myMatrix[1,1]), myMatrix[4,3] - (myMatrix[1,3] * m / myMatrix[1,1]), myMatrix[4,4] - (myMatrix[1,4] * m / myMatrix[1,1])), nrow = 4, ncol = 4, byrow = TRUE)
l4 <- m / myMatrix[1,1]
}
else if ((myMatrix[1,1] != 0) & (m != 0))
{
myMatrix <- matrix(c(myMatrix[1,1], myMatrix[1,2], myMatrix[1,3], myMatrix[1,4], myMatrix[2,1], myMatrix[2,2], myMatrix[2,3], myMatrix[2,4], myMatrix[3,], myMatrix[4,1] - (myMatrix[1,1] * m / myMatrix[1,1]), myMatrix[4,2] - (myMatrix[1,2] * m / myMatrix[1,1]), myMatrix[4,3] - (myMatrix[1,3] * m / myMatrix[1,1]), myMatrix[4,4] - (myMatrix[1,4] * m / myMatrix[1,1])), nrow = 4, ncol = 4, byrow = TRUE)
l4 <- m / myMatrix[1,1]
}
else
{
l4 <- 0
}
n <- myMatrix[4,2]
if ( ((myMatrix[2,2] > 0)&(n > 0)) || ((myMatrix[2,2] < 0)&(n < 0)) )
{
myMatrix <- matrix(c(myMatrix[1,1], myMatrix[1,2], myMatrix[1,3], myMatrix[1,4], myMatrix[2,1], myMatrix[2,2], myMatrix[2,3], myMatrix[2,4], myMatrix[3,], myMatrix[4,1], myMatrix[4,2] - (myMatrix[2,2] * n / myMatrix[2,2]), myMatrix[4,3] - (myMatrix[2,3] * n / myMatrix[2,2]), myMatrix[4,4] - (myMatrix[2,4] * n / myMatrix[2,2])), nrow = 4, ncol = 4, byrow = TRUE)
l5 <- n / myMatrix[2,2]
}
else if ((myMatrix[2,2] != 0) & (n != 0))
{
myMatrix <- matrix(c(myMatrix[1,1], myMatrix[1,2], myMatrix[1,3], myMatrix[1,4], myMatrix[2,1], myMatrix[2,2], myMatrix[2,3], myMatrix[2,4], myMatrix[3,], myMatrix[4,1], myMatrix[4,2] - (myMatrix[2,2] * n / myMatrix[2,2]), myMatrix[4,3] - (myMatrix[2,3] * n / myMatrix[2,2]), myMatrix[4,4] - (myMatrix[2,4] * n / myMatrix[2,2])), nrow = 4, ncol = 4, byrow = TRUE)
l5 <- n / myMatrix[2,2]
}
else
{
l5 <- 0
if ((myMatrix[1,2] == 0) & (myMatrix[2,2] == 0) & (myMatrix[3,2] == 0))
{
l6 <- 0
skipRest <- TRUE
}
}
if (skipRest == FALSE)
{
o <- myMatrix[4,3]
if ( ((myMatrix[3,3] > 0)&(o > 0)) || ((myMatrix[3,3] < 0)&(o < 0)) )
{
myMatrix <- matrix(c(myMatrix[1,1], myMatrix[1,2], myMatrix[1,3], myMatrix[1,4], myMatrix[2,1], myMatrix[2,2], myMatrix[2,3], myMatrix[2,4], myMatrix[3,], myMatrix[4,1], myMatrix[4,2], myMatrix[4,3] - (myMatrix[3,3] * o / myMatrix[3,3]), myMatrix[4,4] - (myMatrix[3,4] * o / myMatrix[3,3])), nrow = 4, ncol = 4, byrow = TRUE)
l6 <- o / myMatrix[3,3]
}
else if ( (myMatrix[3,3] != 0) & (o != 0) )
{
myMatrix <- matrix(c(myMatrix[1,1], myMatrix[1,2], myMatrix[1,3], myMatrix[1,4], myMatrix[2,1], myMatrix[2,2], myMatrix[2,3], myMatrix[2,4], myMatrix[3,], myMatrix[4,1], myMatrix[4,2], myMatrix[4,3] - (myMatrix[3,3] * o / myMatrix[3,3]), myMatrix[4,4] - (myMatrix[3,4] * o / myMatrix[3,3])), nrow = 4, ncol = 4, byrow = TRUE)
l6 <- o / myMatrix[3,3]
}
else
{
l6 <- 0
}
}
print("Your original matrix")
print(saveForLater)
print("Your matrix L is")
lMatrix <- matrix(c(1,0,0,0,l1,1,0,0,l2,l3,1,0,l4,l5,l6,1), nrow = 4, ncol = 4, byrow = TRUE)
print(lMatrix)
print("Your matrix U is")
print(myMatrix)
if (useP == TRUE)
{
print("Your matrix P is")
print(p)
}
}
else
{
print("Cannot factorize. Next time, please enter a matrix that is 3 x 3 or 4 x 4.")
}
}
else
{
print("Cannot factorize. You did not enter a square matrix.")
}
}Now, we will test the function in many ways to show it works.
threeMatrix <- matrix(c(19,2,3,4,2,3,3,1,4), nrow = 3, ncol = 3, byrow = TRUE)
lufactorization(threeMatrix)## [1] "Your original matrix"
## [,1] [,2] [,3]
## [1,] 19 2 3
## [2,] 4 2 3
## [3,] 3 1 4
## [1] "Your matrix L is"
## [,1] [,2] [,3]
## [1,] 1.0000000 0.0000000 0
## [2,] 0.2105263 1.0000000 0
## [3,] 0.1578947 0.4333333 0
## [1] "Your matrix U is"
## [,1] [,2] [,3]
## [1,] 19 2.000000 3.000000
## [2,] 0 1.578947 2.368421
## [3,] 0 0.000000 2.500000threeMatrix <- matrix(c(19,-2,3,-4,2,3,-3,1,4), nrow = 3, ncol = 3, byrow = TRUE)
lufactorization(threeMatrix)## [1] "Your original matrix"
## [,1] [,2] [,3]
## [1,] 19 -2 3
## [2,] -4 2 3
## [3,] -3 1 4
## [1] "Your matrix L is"
## [,1] [,2] [,3]
## [1,] 1.0000000 0.0000000 0
## [2,] -0.2105263 1.0000000 0
## [3,] -0.1578947 0.4333333 0
## [1] "Your matrix U is"
## [,1] [,2] [,3]
## [1,] 19 -2.000000 3.000000
## [2,] 0 1.578947 3.631579
## [3,] 0 0.000000 2.900000threeMatrix <- matrix(c(0,4,1,8,-2,5,3,0,7), nrow = 3, ncol = 3, byrow = TRUE)
lufactorization(threeMatrix)## [1] "Your original matrix"
## [,1] [,2] [,3]
## [1,] 0 4 1
## [2,] 8 -2 5
## [3,] 3 0 7
## [1] "Your matrix L is"
## [,1] [,2] [,3]
## [1,] 1.000 0.0000 0
## [2,] 0.000 1.0000 0
## [3,] 0.375 0.1875 0
## [1] "Your matrix U is"
## [,1] [,2] [,3]
## [1,] 8 -2 5.0000
## [2,] 0 4 1.0000
## [3,] 0 0 4.9375
## [1] "Your matrix P is"
## [,1] [,2] [,3]
## [1,] 0 1 0
## [2,] 1 0 0
## [3,] 0 0 1threeMatrix <- matrix(c(0,4,2,0, 12, 0, 0,0,8), nrow = 3, ncol = 3, byrow = TRUE)
lufactorization(threeMatrix)## [1] " "
## [1] "Your original matrix"
## [,1] [,2] [,3]
## [1,] 0 4 2
## [2,] 0 12 0
## [3,] 0 0 8
## [1] "Your matrix L is"
## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] 0 0 1
## [1] "Your matrix U is"
## [,1] [,2] [,3]
## [1,] 0 4 2
## [2,] 0 12 0
## [3,] 0 0 8
## [1] "This cannot be decomposed any further."## NULLincorrectSize <- matrix(c(1,2,3,3,2,1), nrow = 2, ncol = 3, byrow = TRUE)
lufactorization(incorrectSize)## [1] "Cannot factorize. You did not enter a square matrix."incorrectSize <- matrix(1:49, nrow = 7, ncol = 7, byrow = TRUE)
lufactorization(incorrectSize)## [1] "Cannot factorize. Next time, please enter a matrix that is 3 x 3 or 4 x 4."fourMatrix <- matrix(c(5, 2, 3, 4, 7, 6, 2, 2, 1, 4, 2, 9, 7, 7, 8, 10), nrow = 4, ncol = 4, byrow = TRUE)
lufactorization(fourMatrix)## [1] "Your original matrix"
## [,1] [,2] [,3] [,4]
## [1,] 5 2 3 4
## [2,] 7 6 2 2
## [3,] 1 4 2 9
## [4,] 7 7 8 10
## [1] "Your matrix L is"
## [,1] [,2] [,3] [,4]
## [1,] 1.0 0.0000 0.000000 0
## [2,] 1.4 1.0000 0.000000 0
## [3,] 0.2 1.1250 1.000000 0
## [4,] 1.4 1.3125 1.725806 1
## [1] "Your matrix U is"
## [,1] [,2] [,3] [,4]
## [1,] 5 2.0 3.000 4.00000
## [2,] 0 3.2 -2.200 -3.60000
## [3,] 0 0.0 3.875 12.25000
## [4,] 0 0.0 0.000 -12.01613fourMatrix <- matrix(c(5, 0, 3, 4, 7, 0, 2, 2, 1, 0, 2, 9, 7, 0, 8, 10), nrow = 4, ncol = 4, byrow = TRUE)
lufactorization(fourMatrix)## [1] "Your original matrix"
## [,1] [,2] [,3] [,4]
## [1,] 5 0 3 4
## [2,] 7 0 2 2
## [3,] 1 0 2 9
## [4,] 7 0 8 10
## [1] "Your matrix L is"
## [,1] [,2] [,3] [,4]
## [1,] 1.0 0 0 0
## [2,] 1.4 1 0 0
## [3,] 0.2 0 1 0
## [4,] 1.4 0 0 1
## [1] "Your matrix U is"
## [,1] [,2] [,3] [,4]
## [1,] 5 0 3.0 4.0
## [2,] 0 0 -2.2 -3.6
## [3,] 0 0 1.4 8.2
## [4,] 0 0 3.8 4.4fourMatrix <- matrix(c(-4, -5, 8, 2, 0, -4, 8, 0, 3, 2, -1, -1, 3, 0, 5, 2), nrow = 4, ncol = 4, byrow = TRUE)
lufactorization(fourMatrix)## [1] "Your original matrix"
## [,1] [,2] [,3] [,4]
## [1,] -4 -5 8 2
## [2,] 0 -4 8 0
## [3,] 3 2 -1 -1
## [4,] 3 0 5 2
## [1] "Your matrix L is"
## [,1] [,2] [,3] [,4]
## [1,] 1.000000 0.0000000 0.0 0
## [2,] 0.000000 1.0000000 0.0 0
## [3,] 1.000000 1.7500000 1.0 0
## [4,] -1.333333 -0.6666667 -1.8 1
## [1] "Your matrix U is"
## [,1] [,2] [,3] [,4]
## [1,] 3 2 -1 -1.000000
## [2,] 0 -4 8 0.000000
## [3,] 0 0 -5 3.000000
## [4,] 0 0 0 6.066667
## [1] "Your matrix P is"
## [,1] [,2] [,3] [,4]
## [1,] 0 0 1 0
## [2,] 0 1 0 0
## [3,] 1 0 0 0
## [4,] 0 0 0 1fourMatrix <- matrix(c(0,3,2,1,0,9,3,2,0,0,2,1,0,9,0,3), nrow = 4, ncol = 4, byrow = TRUE)
lufactorization(fourMatrix)## [1] "Your original matrix"
## [,1] [,2] [,3] [,4]
## [1,] 0 3 2 1
## [2,] 0 9 3 2
## [3,] 0 0 2 1
## [4,] 0 9 0 3
## [1] "Your matrix L is"
## [,1] [,2] [,3] [,4]
## [1,] 1 0 0 0
## [2,] 0 1 0 0
## [3,] 0 0 1 0
## [4,] 0 0 0 1
## [1] "Your matrix U is"
## [,1] [,2] [,3] [,4]
## [1,] 0 3 2 1
## [2,] 0 9 3 2
## [3,] 0 0 2 1
## [4,] 0 9 0 3
## [1] "This cannot be decomposed any further."## NULL