(* Evaluate each cell by pressing shift-return or the enter key on the numeric keypad *)

(* Insert matrices by right clicking and selecting insert table/matrix and then selecting matrix *)

a=(({{1, 2}, {2, 4}}));

(* You can also enter matrices as a list of lists with just the keyboard as follows *)

b={{1,2},{3,4}};

(* We can display b as a matrix by typing //MatrixForm *)

b//MatrixForm

(* Bigger matrices are fine too *)

c=(({{1, 2, 3}, {0, 0, 0}, {0, 0, 0}}));

(* You can enter vectors the same way *)

v=(({{1}, {2}}));

(* Vectors and matrices can have variables *)

d=(({{q, w}, {r, t}}));

(Subscript[v, 2])=(({{x}, {y}, {z}}));

(* Find the norm of a vector with norm *)

Norm[v]

(* The norm formula is slightly more complicated if your scalars can be complex instead of real, but that won't affect you for 568 *)

Norm[(Subscript[v, 2])]

(* Take the projection of v=(x,y)
onto the line with direction vector ((Subscript[D, 1]),(Subscript[D, 2])) *)

Projection[{x,y},{(Subscript[D, 1]),(Subscript[D, 2])}, Dot]//MatrixForm

(* Multiply matrices using . *)

a.v //MatrixForm

d.v //MatrixForm

d.d //MatrixForm

(* This is the same as raising d to the second power, which we can do as follows: *)

MatrixPower[d,2]//MatrixForm

(* If there is a general form for a matrix to a power, can solve for that as well.  Here we raise a matrix to the nth power *)

MatrixPower[(({{1, 1}, {0, 1}})),n]//MatrixForm

(* The inverse is the same as taking the matrix to the -1 power, but the shorter way is to just use Inverse *)

MatrixPower[d,-1]//MatrixForm

Inverse[d]//MatrixForm

(* Many matrices are not invertible, but every matrix has a pseudo inverse.  You use the pseudo inverse when finding least squares solutions. *)

PseudoInverse[(({{1, 2}, {1, 2.1}, {2, 3}}))]//MatrixForm

(* Take the transpose of a matrix with Transpose *)

Transpose[d]//MatrixForm

(* Recall:  A matrix is symmetric if it equals its own transpose, so let's test to see if a and b are symmetric *)

a==Transpose[a]

b==Transpose[b]

(* That is, a is symmetric and b is not.  Note we use two = to test for equality.  One equal sign is an assignment *)

(* We can row reduce matrices or find their nullspace *)

(* NullSpace returns a basis for the null space, but transposed from the way you are used to, so we transpose it *)

Transpose[NullSpace[c]]//MatrixForm

(* That is, the vectors (-3;0;1) and (-2;1;0) form a basis for the nullspace of c *)

(* Since b is invertible, the only vector in the null space is the zero vector...  the nullspace is zero dimensional (just one point), so the basis set is empty *)

NullSpace[b]//MatrixForm

(* Row Reducing *)

RowReduce[a]//MatrixForm

RowReduce[b]//MatrixForm

(* Recall we can use row reducing to solve an augmented matrix.  For example: *)

RowReduce[(({{1, 2, 1}, {4, 5, 2}}))]//MatrixForm

(* That is, we have solved the linear system of equations:
 x + 2y == 1 and
4x + 5 y == 2
for x and y, and found x=-(1/3) and y=(2/3) *)

(* Alternatively, we could have just used the solve command *)

Solve[{x + 2y == 1, 4x + 5 y == 2}, {x,y}]

(* Again, symbolic expressions are fine... suppose the right hand sides were q and w *)

Solve[{x + 2y == q, 4x + 5 y == w}, {x,y}]

(* We can find the eigenvalues and eigenvectors of matrices *)

(* There's multiple ways to do this *)

(* First:  From the definition.  We recall the eigenvalues of b are the roots of Det[b-λ IdentityMatrix[2]] *)

Det[b-λ IdentityMatrix[2]]

(* Alternatively, this polynomial is known as the characteristic polynomial of the matrix *)

CharacteristicPolynomial[b, λ]

(* We can then find when the characteristic polynomial is zero *)

(* Some polynomials factor nicely. *)

Factor[(λ^2)-3λ+2]

(* Ours, however, does not *)

Factor[-2-5 λ+(λ^2)]

(* That's okay though.  We can solve for when the characteristic polynomial equals zero. *)

Solve[-2-5 λ+(λ^2)==0, λ]

(* Thus the eigenvalues of b are (1/2) (5±(Sqrt[33])) *)

(* If you only want to know the eigenvalues, just use: *)

Eigenvalues[b]

(* To find just the eigenvectors, use:  (again, due to the way mathematica works, you probably want to take the transpose.)  Here the eigenvectors are (({{-(4/3)+(1/6) (5+(Sqrt[33]))}, {1}})) and (({{-(4/3)+(1/6) (5-(Sqrt[33]))}, {1}})) *)

Transpose[Eigenvectors[b]]//MatrixForm

(* We can check by multiplying the matrix by the vector and comparing it to the eigenvalue times the vector *)

b . (({{-(4/3)+(1/6) (5+(Sqrt[33]))}, {1}})) == (1/2) (5+(Sqrt[33]))  (({{-(4/3)+(1/6) (5+(Sqrt[33]))}, {1}}))

(* Alternatively, use Eigensystem to find both eigenvalues and eigenvectors at the same time *)

{evals, evecs} = Eigensystem[b];

(* The 2nd eigenvalue is: *)

evals[[2]]

(* and it has corresponding eigenvector *)

evecs[[2]]//MatrixForm

(* Remember the determinant of a matrix is the product of the eigenvalues *)

Det[b]

evals[[1]]*evals[[2]]

(* They don't look the same, do they?  We're going to have to simplify. *)

FullSimplify[evals[[1]]*evals[[2]]]

(* This is the same as the determinant *)

(* The trace of b, the sum of the diagonal entries, is (Subscript[b, 11])+(Subscript[b, 22]) *)

Tr[b]

b[[1,1]]+b[[2,2]]

(* The trace is also the sum of the eigenvalues *)

FullSimplify[evals[[1]]+evals[[2]]]

(* To do Gram-Schmidt, use orthogonalize *)

Orthogonalize[{{1,2,3,0},{4,2,1,0},{3,1,2,1}}]

(* That is, the vectors ((1/(Sqrt[14])),(Sqrt[(2/7)]),(3/(Sqrt[14])),0), ((45/(Sqrt[2422])),3 (Sqrt[(2/1211)]),-(19/(Sqrt[2422])),0), etc... are orthogonal vectors, i.e. the dotproduct of any two is zero, with the same span as the original vectors (1,2,3,0), (4,2,1,0) and (3,1,2,1). *)

(* QR decomposition of b gives a matrix Q with orthogonal columns and an upper triangular matrix R such that Q.R=b.  As before, we need to transpose the Q mathematica gives. *)

{Q,R}=QRDecomposition[b];   Q=Transpose[Q];

Q//MatrixForm

R//MatrixForm

Q.R//MatrixForm

(* Q is an orthogonal matrix if and only if Transpose[Q]==Inverse[Q].
Let's check that Q is orthogonal. *)

Transpose[Q]==Inverse[Q]

(* Recall that the determinant of an orthogonal matrix is
±1.  Be careful: the converse is false. *)

Det[Q]

(* For any nonsingular matrix A, there exists matrices P, L and U such that P.A=L.U,
where P is a permutation matrix (think row interchanges),
L is lower triangular with ones on the diagonal and U is upper triangular.  Mathematica returns the lower part of L in the same matrix as U. *)

{mLU, pivots, conditionNum}=LUDecomposition[b]; mLU=Transpose[mLU];

mLU//MatrixForm

(* That is, L=(({{1, 0}, {2, 1}})) and U=(({{1, 3}, {0, -2}})).  Let's check. *)

(({{1, 0}, {2, 1}})).(({{1, 3}, {0, -2}})) //MatrixForm

(* So L.U is b, up to permutation of rows. *)

(* The singular value decomposition decomposes an m x n matrix A into
A=U.D.Transpose[V],
where U and V are orthogonal (unitary in the complex case) and
d is diagonal.  The diagonal entries of d are called the singular values of A. *)

{u,d,v}=SingularValueDecomposition[b];

FullSimplify[u]//MatrixForm

FullSimplify[d]//MatrixForm

FullSimplify[v]//MatrixForm

(* The product u.d.Transpose[v] is the original matrix *)

FullSimplify[u.d.Transpose[v]]//MatrixForm