ASSIGNMENT 1
CUNY 605 FUNDAMENTALS OF COMPUTATIONAL MATHEMATICS - FALL 2017
You can think of vectors representing many dimensions of related information. For instance, Netflix might store all the ratings a user gives to movies in a vector. This is clearly a vector of very large dimensions (in the millions) and very sparse as the user might have rated only a few movies. Similarly, Amazon might store the items purchased by a user in a vector, with each slot or dimension representing a unique product and the value of the slot, the number of such items the user bought. One task that is frequently done in these settings is to find similarities between users. And, we can use dot-product between vectors to do just that. As you know, the dot-product is proportional to the length of two vectors and to the angle between them. In fact, the dot-product between two vectors, normalized by their lengths is called as the cosine distance and is frequently used in recommendation engines.
\[u = \left[\begin{array}{r} 0.5\\ 0.5\\ \end{array} \right] \]
\[v = \left[\begin{array}{r} 3\\ -4\\ \end{array} \right] \]
The dot product of two vectors (if calculated by hand) is: 0.53 + 0.5-4 = -0.5. Let’s confirm this with R.
# Reference: https://stackoverflow.com/questions/4090169/elegant-way-to-check-for-missing-packages-and-install-them
if (!require("geometry")) install.packages("geometry")
library(geometry)
u = c(0.5, 0.5)
v = c(3, -4)
uv.dot <- dot(u,v,d = NULL)
paste0("Dot Product u.v: ", uv.dot)## [1] "Dot Product u.v: -0.5"The lengths of u and v (or otherwise known as ||u|| and ||v||) are calculated as the \(\sqrt{u_{1}^2 + u_{2}^2 + u_{3}^2 + ... u_{n}^2}\) and \(\sqrt{v_{1}^2 + v_{2}^2 + v_{3}^2 + ... v_{n}^2}\) respectively. Again, if we calculated this by hand:
\(||u|| = \sqrt{0.5^2 + 0.5^2} = 0.707106781\)
\(||v|| = \sqrt{3^2 + (-4)^2} = 5\)
Let’s confirm this with R.
# argument x accepts vectors
vector.length <- function(x){
temp <- 0
for (i in 1:length(x)){
temp <- temp + x[i]^2
}
sqrt(temp)
}
u.length <- vector.length(u)
v.length <- vector.length(v)
paste0("Length of Vector U: ", u.length)## [1] "Length of Vector U: 0.707106781186548"paste0("Length of Vector V: ", v.length)## [1] "Length of Vector V: 5"Vectors can perfom addition and scalar multiplication. In this case, this can be re-written as \(3u + (-2)v\).
paste0("3u - 2v = ")## [1] "3u - 2v = "new.uv <- 3*u + (-2)*v
new.uv## [1] -4.5 9.5By using the law of cosine: \(c^2 = a^2 + b^2 - 2ab*cos(\theta)\), we can start to look for the angle between the two vectors. When we substite the different vectors (for each of the sides of the triange), the formula simplifies to \(u*v = ||u|| * ||v|| * cos*(\theta)\). And with this, we can solve for the angle.
theta <- acos(uv.dot/(u.length * v.length))
# Converting radians to degrees
# https://stackoverflow.com/questions/32370485/r-convert-radians-to-degree-degree-to-radians
if (!require("NISTunits")) install.packages("NISTunits", dependencies = TRUE)
library(NISTunits)
deg <- NISTradianTOdeg(theta)
paste0("Angle in Degrees: ", deg, " degrees.")## [1] "Angle in Degrees: 98.1301165223898 degrees."You can use R-markdown to submit your responses to this problem set. If you decide to do it in paper, then please either scan it or take a picture using a smartphone and attach that picture. Please make sure that the picture is legible before submitting.
Set up a system of equations with 3 variables and 3 constraints and solve for x. Please write a function in R that will take two variables (matrix A & constraint vector b) and solve using elimination. Your function should produce the right answer for the system of equations for any 3-variable, 3-equation system. You don’t have to worry about degenerate cases and can safely assume that the function will only be tested with a system of equations that has a solution. Please note that you do have to worry about zero pivots, though. Please note that you should not use the built-in function solve to solve this system or use matrix inverses. The approach that you should employ is to construct an Upper Triangular Matrix and then back-substitute to get the solution. Alternatively, you can augment the matrix A with vector b and jointly apply the Gauss Jordan elimination procedure.
\[ \left[\begin{array}{r} -1.55\\ -0.32\\ 0.95\\ \end{array} \right] \]
\[ \left[\begin{array}{rrr|r} 1 & 1 & 3 & 1 \\ 2 & -1 & 5 & 2 \\ -1 & -2 & 4 & 6 \\ \end{array} \right] \]
We will approach this problem set in 2 ways. We will initially utilize R and approach this question via the Upper Triangular Matrix.
What is the Upper Triangular Matrix?
According to Wikipedia, “a triangular matrix is a special kind of square matrix. A square matrix is called lower triangular if all the entries above the main diagonal are zero. Similarly, a square matrix is called upper triangular if all the entries below the main diagonal are zero…Because matrix equations with triangular matrices are easier to solve, they are very important in numerical analysis.”
We will reference Week 1 Module 2 Reading Assignment for the algorithm that was utilized to help answer this question.
“We can build a simple iterative procedure to systematically convert the coefficient matrix A into the Upper Triangular Form. We will call this procedure Pivot & Multiply.
We will call this function: “solve.matrix(A,b)”, and this function will take argument A (matrix) and b (constraint vector).
# Reference
# http://www.r-tutor.com/r-introduction/matrix
# Defining the matrix and the constraint vector
matrix.A <- matrix(
c(1,1,3,2,-1,5,-1,-2,4), # the data elements
nrow=3, # number of rows
ncol=3, # number of columns
byrow = TRUE) # fill matrix by rows
vector.b <- matrix(
c(1,2,6), # the data elements
nrow=3, # number of rows
ncol=1, # number of columns
byrow = TRUE) # fill matrix by rows
# Argument A = matrix
# Argument b = constraint vectorLet’s define the function:
solve.matrix <- function(A,b){
# This function only works via the Upper Triangular Matrix
# The upper triangular matrix is if all the entries below the main diagonal are zero.
# Argument A is the matrix
# Combine A with b to make augmented matrix
Aug.A <- cbind(A,b)
# The question stems asks that a function be specifically be built for a 3 x 3 matrix. We will assume that the argument A is a 3 x 3 matrix.
# This means that we will need to zero out [2,1], [2,3], and [3,2]
# Start with row 1 of the co-efficient matrix
# Pivot: Evaluate the first non-zero element
pivot <- Aug.A[1,1]
# Multiplier: The element being eliminated divided by the Pivot
# Will need to zero out [2,1] by dividing this number by the pivot.
multiplier <- Aug.A[2,1]/pivot
# Subtract Multiplier times row n from row n+1 like the Gauss Jordan elimination
Aug.A[2, ] <- Aug.A[2,] - Aug.A[1,]*multiplier
# Again, will repeat the previous steps for [3,1]
multiplier <- Aug.A[3,1]/pivot
Aug.A[3, ] <- Aug.A[3, ] - Aug.A[1, ]*multiplier
# And one last time, will repeat these steps for [3,2]
# Given that we are now in column 2, we will need to find a new pivot
# And since [2,1] already contains a zero, we naturally choose [2,2] as the next pivot point as that is a non-zero number.
pivot <- Aug.A[2,2]
multiplier <- Aug.A[3,2]/pivot
Aug.A[3, ] <- Aug.A[3, ] - Aug.A[2, ]*multiplier
# Now that we have successfully created the upper triangular matrix, we can now solve for each n variable
# We'll go backwards.
x3 <- Aug.A[3,4]/Aug.A[3,3]
x2 <- (Aug.A[2,4] - Aug.A[2,3]*x3)/Aug.A[2,2]
x1 <- (Aug.A[1,4] - Aug.A[1,3]*x3 - Aug.A[1,2]*x2)/Aug.A[1,1]
ans.vector <- c(x1,x2,x3)
ans.vector
}Let’s confirm this with a back-substituion:
ans.vec <- solve.matrix(matrix.A, vector.b)
x1 <- ans.vec[1]
x2 <- ans.vec[2]
x3 <- ans.vec[3]
ans1 <- x1 + x2 + 3*x3
ans2 <- 2*x1 -x2 + 5*x3
ans3 <- -1*x1 - 2*x2 + 4*x3
paste0("x1 + x2 + 3*x3 = ", ans1)## [1] "x1 + x2 + 3*x3 = 1"paste0("2*x1 - x2 + 5*x3 = ", ans2)## [1] "2*x1 - x2 + 5*x3 = 2"paste0("-1*x1 - 2*x2 + 4*x3 = ", ans3)## [1] "-1*x1 - 2*x2 + 4*x3 = 6"Now that we had confirmed our answers via the Upper Triangular Matrix method, we will now proceed by hand to answer the same question, but instead with the Gauss Jordan Elimination Method.
For a refresher, the Gauss Jordan Elimination method entails a reduction of the matrix to its row echelon form via these pathways:
Let’s start by solving this problem by hand and showing the work:
\[ \left[\begin{array}{rrr|r} 1 & 1 & 3 & 1 \\ 0 & -3 & -1 & 0 \\ -1 & -2 & 4 & 6 \\ \end{array} \right] \]
\[ \left[\begin{array}{rrr|r} 1 & 1 & 3 & 1 \\ 0 & -3 & -1 & 0 \\ 0 & -1 & 7 & 7 \\ \end{array} \right] \]
\[ \left[\begin{array}{rrr|r} 1 & 0 & 10 & 8 \\ 0 & -3 & -1 & 0 \\ 0 & -1 & 7 & 7 \\ \end{array} \right] \]
\[ \left[\begin{array}{rrr|r} 1 & 0 & 10 & 8 \\ 0 & 0 & -22 & -21 \\ 0 & -1 & 7 & 7 \\ \end{array} \right] \]
\[ \left[\begin{array}{rrr|r} 1 & 0 & 10 & 8 \\ 0 & 0 & -22 & -21 \\ 0 & 1 & -7 & -7 \\ \end{array} \right] \]
\[ \left[\begin{array}{rrr|r} 1 & 0 & 10 & 8 \\ 0 & 1 & -7 & -7 \\ 0 & 0 & -22 & -21 \\ \end{array} \right] \]
\[ \left[\begin{array}{rrr|r} 1 & 0 & 10 & 8 \\ 0 & 1 & -7 & -7 \\ 0 & 0 & 1 & 21/22 \\ \end{array} \right] \]
\[ \left[\begin{array}{rrr|r} 1 & 0 & 10 & 8 \\ 0 & 1 & 0 & -7/22 \\ 0 & 0 & 1 & 21/22 \\ \end{array} \right] \]
\[ \left[\begin{array}{rrr|r} 1 & 0 & 0 & -17/11 \\ 0 & 1 & 0 & -7/22 \\ 0 & 0 & 1 & 21/22 \\ \end{array} \right] \]
The above is the solution via the Gauss Jordan elimination. We will now back-substitute to get the solution.
x1 <- -17/11
x2 <- -7/22
x3 <- 21/22
ans1 <- x1 + x2 + 3*x3
ans2 <- 2*x1 -x2 + 5*x3
ans3 <- -1*x1 - 2*x2 + 4*x3
paste0("x1 + x2 + 3*x3 = ", ans1)## [1] "x1 + x2 + 3*x3 = 1"paste0("2*x1 - x2 + 5*x3 = ", ans2)## [1] "2*x1 - x2 + 5*x3 = 2"paste0("-1*x1 - 2*x2 + 4*x3 = ", ans3)## [1] "-1*x1 - 2*x2 + 4*x3 = 6"As you can see, we again confirm this answer.