딥스탯 2018. 10. 4. 15:22
11_Basic_linear_algebra

Basic linear algebra in Julia

Author: Andreas Noack Jensen (MIT) (http://www.econ.ku.dk/phdstudent/noack/) (with edits from Jane Herriman)

Reference

https://github.com/JuliaComputing/JuliaBoxTutorials/tree/master/introductory-tutorials/intro-to-julia (github : JuliaComputing/JuliaBoxTutorials/introductory-tutorials/intro-to-julia/)

Topics:

  1. Basic operations
  2. The LinearAlgebra library
  3. Exercises

Series

First let's define a random matrix

In [1]:
A = rand(1:4,4,4)
Out[1]:
4×4 Array{Int64,2}:
 4  4  1  1
 1  4  3  3
 4  4  4  4
 3  4  4  2

Define a vector of ones

In [2]:
x = fill(1.0, (4,)) # = fill(1.0, 4)
Out[2]:
4-element Array{Float64,1}:
 1.0
 1.0
 1.0
 1.0

Notice that A has type Array{Int64,2} but x has type Array{Float64,1}. Julia defines the aliases Vector{Type}=Array{Type,1} and Matrix{Type}=Array{Type,2}.

Many of the basic operations are the same as in other languages

Multiplication

In [3]:
b = A*x
Out[3]:
4-element Array{Float64,1}:
 10.0
 11.0
 16.0
 13.0

Transposition

As in other languages A' is the conjugate transpose, or adjoint

In [4]:
A'
Out[4]:
4×4 LinearAlgebra.Adjoint{Int64,Array{Int64,2}}:
 4  1  4  3
 4  4  4  4
 1  3  4  4
 1  3  4  2

and we can get the transpose with

In [5]:
transpose(A)
Out[5]:
4×4 LinearAlgebra.Transpose{Int64,Array{Int64,2}}:
 4  1  4  3
 4  4  4  4
 1  3  4  4
 1  3  4  2

Transposed multiplication

Julia allows us to write this without *

In [6]:
A'A
Out[6]:
4×4 Array{Int64,2}:
 42  48  35  29
 48  64  48  40
 35  48  42  34
 29  40  34  30

Solving linear systems

The problem Ax = b for square is solved by the \ function.

In [7]:
A\b
Out[7]:
4-element Array{Float64,1}:
 1.0
 1.0
 1.0
 1.0
In [8]:
A^-1 * b
Out[8]:
4-element Array{Float64,1}:
 0.9999999999999996
 1.0000000000000004
 1.0               
 0.9999999999999991

A\b gives us the least squares solution if we have an overdetermined linear system (a "tall" matrix)

In [9]:
Atall = rand(4, 2)
Out[9]:
4×2 Array{Float64,2}:
 0.620125  0.964271
 0.435129  0.392623
 0.352525  0.846371
 0.888477  0.154265
In [10]:
Atall\b
Out[10]:
2-element Array{Float64,1}:
 12.558191782805045
  7.886509151359657
In [11]:
(Atall'Atall)^-1 * Atall'b
Out[11]:
2-element Array{Float64,1}:
 12.55819178280504 
  7.886509151359654

and the minimum norm least squares solution if we have a rank-deficient least squares problem

In [12]:
v1 = rand(4)
v2 = rand(4)
rankdef = hcat(v1, v1, v2)
Out[12]:
4×3 Array{Float64,2}:
 0.791797  0.791797  0.631865
 0.132592  0.132592  0.468963
 0.421476  0.421476  0.177085
 0.967094  0.967094  0.202275
In [13]:
rankdef\b
Out[13]:
3-element Array{Float64,1}:
  5.49389774957956 
  5.493897749579558
 11.26631347260693 

Julia also gives us the minimum norm solution when we have an underdetermined system (a "short" matrix)

In [14]:
bshort = rand(2)
Ashort = rand(2, 3)
Out[14]:
2×3 Array{Float64,2}:
 0.913861   0.382365  0.542242
 0.0607539  0.420378  0.750087
In [15]:
Ashort\bshort
Out[15]:
3-element Array{Float64,1}:
 0.279885521561581  
 0.37741831223794353
 0.6366857785894317 

The LinearAlgebra library

While much of linear algebra is available in Julia by default (as shown above), there's a standard library named LinearAlgebra that brings in many more relevant names and functions. In particular, it provides factorizations and some structured matrix types. As with all packages, you can bring these additional features into your session with a using LinearAlgebra.

Exercises

10.1

Take the inner product (or "dot" product) of a vector v with itself and assign it to variable dot_v.

In [16]:
v = [1,2,3]
Out[16]:
3-element Array{Int64,1}:
 1
 2
 3
In [17]:
dot_v = v'v
Out[17]:
14
In [18]:
using LinearAlgebra
In [19]:
dot_v = dot(v,v)
Out[19]:
14
In [20]:
dot_v = v⋅v # \cdot<tab>
Out[20]:
14
In [21]:
dot_v = (v,v) # \cdot<tab>
Out[21]:
14
In [22]:
@assert dot_v == 14

10.2

Take the outer product of a vector v with itself and assign it to variable cross_v

In [23]:
using LinearAlgebra
In [24]:
cross_v = cross(v,v)
Out[24]:
3-element Array{Int64,1}:
 0
 0
 0
In [25]:
cross_v = v×v # \times<tab>
Out[25]:
3-element Array{Int64,1}:
 0
 0
 0
In [26]:
cross_v = ×(v,v) # \times<tab>
Out[26]:
3-element Array{Int64,1}:
 0
 0
 0
In [27]:
@assert cross_v == [0, 0, 0]