Basic Matrix Algebra

This notebook presents some basic linear algebra in Julia.


In [1]:
using Dates, Printf, LinearAlgebra

include("printmat.jl")


Out[1]:
printwhere (generic function with 1 method)

Adding and Multiplying: A Matrix and a Scalar

With a matrix $A$ and a scalar $c$, do

  1. A*c (textbook: $Ac$) to multiply each element of $A$ by $c$

  2. A .+ c (textbook: $A+cJ$, where $J$ is a matrix of ones) to add $c$ to each element of $A$, and similarly A .- c ($A-c$)

Watch out when the number comes first: 2.+A is not allowed since it is ambiguous. However, 2.0.+A and 2 .+ A both work.


In [2]:
A = [1 3;3 4]
c = 10

println("A:")
printmat(A)
println("c:")
printmat(c)

println("A*c:")
printmat(A*c)

println("A .+ c:")
printmat(A .+ c)          #notice the dot in .+


A:
         1         3
         3         4

c:
        10

A*c:
        10        30
        30        40

A .+ c:
        11        13
        13        14

Adding and Multiplying Two Matrices

With two matrices of the same dimensions ($A$ and $B$), do

A+B (textbook: $A+B$) to add them (element by element), and similarly A-B (textbook: $A-B$).

Multiplying matrices ($A$ and $B$) of conformable dimensions

A*B (textbook: $AB$)


In [3]:
A = [1 3;3 4]               #A and B are 2x2 matrices
B = [1 2;3 -2]
println("A:")
printmat(A)
println("B:")
printmat(B)

println("A+B:")
printmat(A+B)

println("A*B:")
printmat(A*B)


A:
         1         3
         3         4

B:
         1         2
         3        -2

A+B:
         2         5
         6         2

A*B:
        10        -4
        15        -2

Transpose

You can transpose a numerical matrix A by A'.

Notice that (in Julia) A and B=A' share the same elements (changing one changes the other). If you want an independent copy, use B=copy(A').

For an array of other elements (for instance, strings), use permutedims(A) to swap the dimensions.


In [4]:
A = [1 2 3;4 5 6]
println("A: ")
printmat(A)
println("A': ")
printmat(A')


A: 
         1         2         3
         4         5         6

A': 
         1         4
         2         5
         3         6

Matrix Inverse

A matrix inverse of an $nxn$ matrix $A$:

inv(A) or A^(-1) (textbook: $A^{-1}$)

The inverse is such that $AA^{-1}=I$ and $A^{-1}A=I$


In [5]:
A = [1 3;3 4]
println("A:")
printmat(A)

println("inv(A):")
printmat(inv(A))

println("inv(A)*A:")
printmat(inv(A)*A)


A:
         1         3
         3         4

inv(A):
    -0.800     0.600
     0.600    -0.200

inv(A)*A:
     1.000    -0.000
     0.000     1.000

Vectors: Inner and Outer Products

There are several different ways to think about a vector in mathematics: as a $K \times 1$ matrix (a column vector), a $1 \times K$ matrix (a row vector) or just a flat $K-$vector. Julia uses flat vectors but they are mostly interchangable with column vectors.

The inner product of two (column) vectors with $k$ elements is calculated as x'z or dot(x,y) (textbook: $x'z$ or $x \cdot z$) to get a scalar. (You can also use or x⋅y where the dot is obtained by \cdot + TAB, but that is sometimes hard to distinguish from or x.y)

In contrast, the outer of two (column) vectors with $k$ elements is calculated (textbook: $xz'$) to get a $k\times k$ matrix.


In [6]:
x = [10,11]                  #[10;11] gives the same
z = [2,5]
println("x and z")
printmat([x z])

println("x'z: ")
printlnPs(x'z)               #dot(x,z) gives the same

println("x*z':")
printmat(x*z')


x and z
        10         2
        11         5

x'z: 
        75
x*z':
        20        50
        22        55

Vectors: Quadratic Forms

A quadratic form ($A$ is an $n \times n$ matrix and x is an $n$ vector): x'A*x (textbook: $x'Ax$) to get a scalar.

(In Julia 1.5+ there is also the form dot(x,A,x))


In [7]:
A = [1 3;3 4]
x = [10,11]

println("x:")
printmat(x)
println("A:")
printmat(A)

println("x'A*x: ")
printlnPs(x'A*x)


x:
        10
        11

A:
         1         3
         3         4

x'A*x: 
      1244

Vectors: Extracting Vectors from Matrices

Notice that A[1,:] and A[:,1] both give flat vectors. In case you want a row vector use A[1:1,:].


In [8]:
A[1,:]


Out[8]:
2-element Array{Int64,1}:
 1
 3

OLS Notation

$X'X$ or $\sum\nolimits_{t=1}^{T}x_{t}x_{t}^{\prime}$?

Let $x_t$ be a (column) vector with values of $K$ regressors for observation $t$. Then $x_{t}x_{t}^{\prime}$ is the outer product (a $K\times K$ matrix) and $\sum\nolimits_{t=1}^{T}x_{t}x_{t}^{\prime}$ is just the sum (of each element) across the $T$ observations.

We can calculate the same thing by (a) letting $X$ be a $T\times K$ matrix with $x_{t}^{\prime}$ in row $t$ and (b) then do $X'X$.


In [9]:
x₁ = [1,-1]     #a (column) vector
x₂ = [1,0]
x₃ = [1,1.0]

X = [x₁';x₂';x₃']

println("X")
printmat(X)

(T,K) = (size(X,1),size(X,2))

Sxx1 = x₁*x₁' + x₂*x₂' + x₃*x₃'      #just to illustrate

Sxx2 = zeros(K,K)
for t = 1:T
    #global Sxx2                   #only needed in REPL/script
    Sxx2 = Sxx2 + X[t,:]*X[t,:]'   #X[t,:] becomes a flat vector
end

Sxx3 = X'X

println("sum of outer products, three versions")
printmat(Sxx1)
printmat(Sxx2)
printmat(Sxx3)


X
     1.000    -1.000
     1.000     0.000
     1.000     1.000

sum of outer products, three versions
     3.000     0.000
     0.000     2.000

     3.000     0.000
     0.000     2.000

     3.000     0.000
     0.000     2.000


In [ ]: