library('matrixcalc')

1.1 and 1.2 \[ \text{Given A is } m \times n, \space m \neq n: \\ dim(A^\intercal A) = (n \times m) \times (m \times n) = n \times n \stackrel{?}{=} m \times m = (m \times n) \times (n \times m) = dim(AA^\intercal) \\ \text{But we are given } m \neq n, \space \therefore \space A^\intercal A \neq AA^\intercal\\ \text{Given A is } m \times m, \space \text{with components } a,b,c,d \in \mathbb{C}: \\ A = \begin{bmatrix} a & b \\ c & d \\ \end{bmatrix} \\ A^\intercal A = \begin{bmatrix} a & c \\ b & d \\ \end{bmatrix} \times \begin{bmatrix} a & b \\ c & d \\ \end{bmatrix} = \begin{bmatrix} a^2 + c^2 & ab + cd \\ ab + cd & b^2 + d^2 \\ \end{bmatrix} \\ AA^\intercal = \begin{bmatrix} a & b \\ c & d \\ \end{bmatrix} \times \begin{bmatrix} a & c \\ b & d \\ \end{bmatrix} = \begin{bmatrix} a^2 + b^2 & ac + bd \\ ac + bd & c^2 + d^2 \\ \end{bmatrix} \\ A^\intercal A = AA^\intercal \text{ if :} \\ \begin{cases} a^2 + c^2 = a^2 + b^2 \implies c^2 = b^2 \\ ab + cd = ac + bd \implies a(b-c) = d(b - c)\\ b^2 + d^2 = c^2 + d^2 \implies b^2 = c^2\\ \end{cases} \\ \therefore A^\intercal A = AA^\intercal \text{ if } b = c \implies A \text{ is symmetric} \] 1.1 and 1.2 examples

u = c(.5, .5)
v = c(3, -4)
w = c(0, -5)
z = c(10, 0)
A = cbind(u,v,z)
B = cbind(w,z)
C = cbind(z,w)
print('A, not square')
## [1] "A, not square"
A
##        u  v  z
## [1,] 0.5  3 10
## [2,] 0.5 -4  0
print('AAt')
## [1] "AAt"
A %*% t(A)
##        [,1]   [,2]
## [1,] 109.25 -11.75
## [2,] -11.75  16.25
print('AtA')
## [1] "AtA"
t(A) %*% A
##      u    v   z
## u  0.5 -0.5   5
## v -0.5 25.0  30
## z  5.0 30.0 100
print('B, square but not symmetric')
## [1] "B, square but not symmetric"
B
##       w  z
## [1,]  0 10
## [2,] -5  0
print('BBt')
## [1] "BBt"
B %*% t(B)
##      [,1] [,2]
## [1,]  100    0
## [2,]    0   25
print('BtB')
## [1] "BtB"
t(B) %*% B
##    w   z
## w 25   0
## z  0 100
print('C, square and symmetric')
## [1] "C, square and symmetric"
C
##       z  w
## [1,] 10  0
## [2,]  0 -5
print('CCt')
## [1] "CCt"
C %*% t(C)
##      [,1] [,2]
## [1,]  100    0
## [2,]    0   25
print('CtC')
## [1] "CtC"
t(C) %*% C
##     z  w
## z 100  0
## w   0 25

2

A = cbind(c(1,2,-1), c(1,-1,-2), c(3,5,4))
b = c(1,2,6)
A2 = cbind(c(0,2,-1), c(1,-1,-2), c(3,5,4))

LU_factor = function(A, verbose=FALSE) {
  # LU_factorization as from https://www.youtube.com/watch?v=UlWcofkUDDU
  # setup L and U
  # Whatever you do to clear out row,col in U, put the negative in row,col in L
  # assume no row swaps needed
  dims_A = dim(A)
  L = diag(dims_A[1])
  U = A
  #row_lim = 1
  
  for (col_ix in 1:ncol(A)) {
    for (row_ix in 1:(nrow(A)-1)) {
      # if you are on the diagonal
      if (row_ix == col_ix) {
        clear_row = U[row_ix,]
        # use clear_row to clear out each nonzero entry below diagonal
        for (other_row_ix in (row_ix+1):nrow(A)) {
          if (U[other_row_ix,col_ix] != 0) {
            clear_multi = U[other_row_ix,col_ix]/U[row_ix, col_ix]
            U[other_row_ix,] = U[other_row_ix,] - clear_row * clear_multi
            L[other_row_ix, col_ix] = clear_multi
            #print(L)
          }
          if (verbose) {
            print(L)
            print(U)
            }
        }
        break
      }
    }
  }
  return(list(L,U))
}

Test some examples

A = cbind(c(1,2,-1), c(1,-1,-2), c(3,5,4))
A
##      [,1] [,2] [,3]
## [1,]    1    1    3
## [2,]    2   -1    5
## [3,]   -1   -2    4
A
##      [,1] [,2] [,3]
## [1,]    1    1    3
## [2,]    2   -1    5
## [3,]   -1   -2    4
LU_list = LU_factor(A, verbose=FALSE)
L = as.matrix(LU_list[[1]])
U = as.matrix(LU_list[[2]])
L
##      [,1]      [,2] [,3]
## [1,]    1 0.0000000    0
## [2,]    2 1.0000000    0
## [3,]   -1 0.3333333    1
U
##      [,1] [,2]      [,3]
## [1,]    1    1  3.000000
## [2,]    0   -3 -1.000000
## [3,]    0    0  7.333333
L %*% U
##      [,1] [,2] [,3]
## [1,]    1    1    3
## [2,]    2   -1    5
## [3,]   -1   -2    4
LU_decomp = lu.decomposition(A)
LU_decomp$L
##      [,1]      [,2] [,3]
## [1,]    1 0.0000000    0
## [2,]    2 1.0000000    0
## [3,]   -1 0.3333333    1
LU_decomp$U
##      [,1] [,2]      [,3]
## [1,]    1    1  3.000000
## [2,]    0   -3 -1.000000
## [3,]    0    0  7.333333
LU_decomp$L %*% LU_decomp$U
##      [,1] [,2] [,3]
## [1,]    1    1    3
## [2,]    2   -1    5
## [3,]   -1   -2    4

Test another example

a<-matrix(c(2,4,-4,1,-4,3,-6,-9,5), byrow = T, nrow=3)
soln<-lu.decomposition(a)
soln[[1]]%*%soln[[2]]
##      [,1] [,2] [,3]
## [1,]    2    4   -4
## [2,]    1   -4    3
## [3,]   -6   -9    5
my_soln = LU_factor(a)
my_soln[[1]]%*%my_soln[[2]]
##      [,1] [,2] [,3]
## [1,]    2    4   -4
## [2,]    1   -4    3
## [3,]   -6   -9    5