In [1]:
## Lecture 4, sparse matrix operations.

#= We are going to look at the Candyland matrix!

See some analysis online:
* <http://www.lscheffer.com/CandyLand.htm>
* <http://datagenetics.com/blog/december12011/index.html>

We are studying a slightly different set of rules
with respect to the sticky states. In one set
of rules (on those pages), a sticky state just holds
your character determinstically for one step. In our
set of rules, a sticky state holds you until you draw
a specific color. So if you see some differences,
that is why.
=#
In [2]:
## We have build the Markov transition matrix for you!
using DelimitedFiles
using SparseArrays
data = readdlm("candyland-matrix.csv",',')
TI = Int.(data[:,1])
TJ = Int.(data[:,2])
TV = data[:,3]
T = sparse(TI,TJ, TV, 140,140)
coords = readdlm("candyland-coords.csv",',')
x = coords[:,1]
y = coords[:,2]
cells = Int.(readdlm("candyland-cells.csv",','))
start,dest = cells[1:2]
Out[2]:
2-element Array{Int64,1}:
 140
 134
In [3]:
##
sum(T,dims=1)
Out[3]:
1×140 Array{Float64,2}:
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  …  0.0  1.0  1.0  1.0  1.0  1.0  1.0
In [4]:
##
findall(vec(sum(T,dims=1)) .< 0.9)
Out[4]:
1-element Array{Int64,1}:
 134
In [5]:
## Let's see how Julia stores thigns
typeof(T)
Out[5]:
SparseMatrixCSC{Float64,Int64}
In [6]:
## What can we do with T?
fieldnames(typeof(T))
(:m, :n, :colptr, :rowval, :nzval)
Out[6]:
(:m, :n, :colptr, :rowval, :nzval)
In [7]:
##
T.colptr
Out[7]:
141-element Array{Int64,1}:
    1
   19
   37
   55
   73
   74
   92
  110
  128
  146
  164
  182
  200
    ⋮
 2210
 2220
 2229
 2237
 2244
 2244
 2245
 2246
 2249
 2252
 2255
 2273
In [8]:
##
T.nzval
Out[8]:
2272-element Array{Float64,1}:
 0.09375 
 0.09375 
 0.09375 
 0.09375 
 0.09375 
 0.09375 
 0.0625  
 0.015625
 0.0625  
 0.0625  
 0.046875
 0.046875
 0.0625  
 ⋮       
 0.0625  
 0.0625  
 0.015625
 0.0625  
 0.0625  
 0.046875
 0.046875
 0.015625
 0.015625
 0.015625
 0.015625
 0.015625
In [9]:
##
T.rowval
Out[9]:
2272-element Array{Int64,1}:
   2
   3
   4
   5
   6
   7
   8
   9
  10
  11
  12
  13
  14
   ⋮
   7
   8
   9
  10
  11
  12
  13
  20
  42
  69
  92
 102
In [10]:
## Why sparse matrices?
# The first reason is memory!
print(varinfo(r"T"))
A = Array(T)
print(varinfo(r"A")) # so it takes about 10% of the memory!

| name |       size | summary                                |
|:---- | ----------:|:-------------------------------------- |
| T    | 36.758 KiB | 140×140 SparseMatrixCSC{Float64,Int64} |
| TI   | 17.789 KiB | 2272-element Array{Int64,1}            |
| TJ   | 17.789 KiB | 2272-element Array{Int64,1}            |
| TV   | 17.789 KiB | 2272-element Array{Float64,1}          |
| name |        size | summary                  |
|:---- | -----------:|:------------------------ |
| A    | 153.164 KiB | 140×140 Array{Float64,2} |
In [11]:
## The second reason is speed as you'll see on the homework.
In [12]:
## The downside is that some operations need a lot more care.
In [13]:
## Let's write our own matvec given the raw arrays
function mymatvec(TI,TJ,TV,x,m)
  y = zeros(m)
  for nzi in 1:length(TI)
    i = TI[nzi]
    j = TJ[nzi]
    v = TV[nzi]
    y[i] += v*x[j]
  end
  # A fancy way
  #for (i,j,v) in zip(TI,TJ,TV)
  #  y[i] += v*x[j]
  # end
  return y
end
x = randn(140)
mymatvec(TI,TJ,TV,x,140) == T*x
Out[13]:
true
In [14]:
##
y = mymatvec(TI,TJ,TV,x,140)
z = T*x
Out[14]:
140-element Array{Float64,1}:
  0.020876579243726497
 -0.15226256478824346 
 -0.1753229896049116  
 -0.2425933171008695  
 -0.31776504704772707 
 -0.31776504704772707 
 -0.12826052337856453 
 -0.2568626891779052  
  0.15028666030108256 
 -0.39037630305389065 
 -0.5015320253900685  
 -0.4528017515027682  
 -0.49790114614004083 
  ⋮                   
  0.6798451202517689  
  0.6606617566383506  
  0.5657793460405882  
  0.34927363571129333 
  0.3196337576595978  
  1.7525247067902565  
 -0.39291721229314464 
  1.583282958194908   
 -1.277960226924249   
  0.052839856405129226
  1.2933476267791673  
  0.0
In [15]:
##
## Let's see how Julia stores thigns
typeof(T)
Out[15]:
SparseMatrixCSC{Float64,Int64}
In [16]:
## What can we do with T?
fieldnames(typeof(T))
Out[16]:
(:m, :n, :colptr, :rowval, :nzval)
In [17]:
##
T.colptr
Out[17]:
141-element Array{Int64,1}:
    1
   19
   37
   55
   73
   74
   92
  110
  128
  146
  164
  182
  200
    ⋮
 2210
 2220
 2229
 2237
 2244
 2244
 2245
 2246
 2249
 2252
 2255
 2273
In [18]:
##
T.nzval
Out[18]:
2272-element Array{Float64,1}:
 0.09375 
 0.09375 
 0.09375 
 0.09375 
 0.09375 
 0.09375 
 0.0625  
 0.015625
 0.0625  
 0.0625  
 0.046875
 0.046875
 0.0625  
 ⋮       
 0.0625  
 0.0625  
 0.015625
 0.0625  
 0.0625  
 0.046875
 0.046875
 0.015625
 0.015625
 0.015625
 0.015625
 0.015625
In [19]:
##
T.rowval
Out[19]:
2272-element Array{Int64,1}:
   2
   3
   4
   5
   6
   7
   8
   9
  10
  11
  12
  13
  14
   ⋮
   7
   8
   9
  10
  11
  12
  13
  20
  42
  69
  92
 102
In [20]:
##
function cscmatvec(colptr,rowval,nzval,x,m)
  y = zeros(m)
  for j=1:length(colptr)-1 # for each column ...
    for nzi=colptr[j]:colptr[j+1]-1 # for each entry in the column
      i = rowval[nzi]
      v = nzval[nzi]
      y[i] += v*x[j]
    end
  end
  return y
end
cscmatvec(T.colptr,T.rowval,T.nzval,x,140) == T*x
Out[20]:
true
In [21]:
## Note that getting a random column in T is fast
# Note, run this a few times, the first time
# in Julia tends to be unreliable.
@time T[:,(rand(1:size(T,2)))]

  0.062905 seconds (94.56 k allocations: 4.546 MiB, 14.43% gc time)
Out[21]:
140-element SparseVector{Float64,Int64} with 14 stored entries:
  [9  ]  =  0.015625
  [20 ]  =  0.015625
  [42 ]  =  0.015625
  [69 ]  =  0.015625
  [92 ]  =  0.015625
  [102]  =  0.015625
  [127]  =  0.09375
  [128]  =  0.09375
  [129]  =  0.09375
  [130]  =  0.09375
  [131]  =  0.09375
  [132]  =  0.09375
  [133]  =  0.0625
  [134]  =  0.28125
In [22]:
## But getting a random row in T is slow
@time T[rand(1:size(T,1)),:]

  0.046951 seconds (90.38 k allocations: 4.309 MiB)
Out[22]:
140-element SparseVector{Float64,Int64} with 11 stored entries:
  [115]  =  0.0625
  [116]  =  0.0625
  [118]  =  0.0625
  [119]  =  0.0625
  [120]  =  0.0625
  [121]  =  0.09375
  [122]  =  0.09375
  [123]  =  0.09375
  [124]  =  0.09375
  [125]  =  0.09375
  [126]  =  0.09375
In [23]:
## Okay, maybe we need a bigger matrix
S = sprandn(100000,100000,10/100000)
Out[23]:
100000×100000 SparseMatrixCSC{Float64,Int64} with 999746 stored entries:
  [416   ,      1]  =  -1.17697
  [14586 ,      1]  =  -0.0714233
  [23673 ,      1]  =  -1.00384
  [49266 ,      1]  =  0.199751
  [53418 ,      1]  =  0.635684
  [67465 ,      1]  =  1.43758
  [69709 ,      1]  =  0.95576
  [77109 ,      1]  =  -0.500347
  [87284 ,      1]  =  -0.654499
  [5767  ,      2]  =  1.45467
  [18826 ,      2]  =  -1.67098
  [19650 ,      2]  =  -0.24733
  ⋮
  [49486 ,  99999]  =  0.0996024
  [58370 ,  99999]  =  -0.210823
  [59800 ,  99999]  =  0.419807
  [60788 ,  99999]  =  0.221586
  [75571 ,  99999]  =  0.183908
  [82751 ,  99999]  =  -1.52789
  [7904  , 100000]  =  -0.175436
  [8529  , 100000]  =  -0.131205
  [17768 , 100000]  =  -1.70744
  [33942 , 100000]  =  -1.29998
  [73226 , 100000]  =  -0.177817
  [75057 , 100000]  =  1.84214
In [24]:
##Note that getting a random column in S is fast
@time S[:,rand(1:size(T,2))]

  0.000011 seconds (8 allocations: 512 bytes)
Out[24]:
100000-element SparseVector{Float64,Int64} with 7 stored entries:
  [9445  ]  =  0.0923899
  [14186 ]  =  0.498531
  [46278 ]  =  -0.591371
  [78435 ]  =  -0.332403
  [90408 ]  =  1.12001
  [95259 ]  =  0.820223
  [96095 ]  =  -0.373651
In [25]:
## But getting a random row in S is slow
@time S[rand(1:size(T,1)),:]

  0.000732 seconds (14 allocations: 928 bytes)
Out[25]:
100000-element SparseVector{Float64,Int64} with 9 stored entries:
  [23991 ]  =  -0.395106
  [33141 ]  =  0.0531024
  [46387 ]  =  0.0121212
  [50005 ]  =  1.64379
  [70507 ]  =  0.0605798
  [82804 ]  =  1.85504
  [92190 ]  =  -0.160045
  [93535 ]  =  -1.57911
  [98611 ]  =  -0.677457
In [26]:
## Also note that there is NO way we could possibly store S with dense
length(S)*8/1e9 # it would take 80GB of memory, which is 10x my Mac.
Out[26]:
80.0
In [27]:
## And we can go bigger
S = sprandn(10000000,10000000,10/10000000)
Out[27]:
10000000×10000000 SparseMatrixCSC{Float64,Int64} with 99996212 stored entries:
  [573262  ,        1]  =  -0.542425
  [1895088 ,        1]  =  -0.461467
  [2120880 ,        1]  =  -0.226537
  [2129636 ,        1]  =  1.43839
  [2343463 ,        1]  =  0.84837
  [2806957 ,        1]  =  -1.79241
  [3110223 ,        1]  =  2.17915
  [3487095 ,        1]  =  0.0806138
  [3644555 ,        1]  =  0.351717
  [4313391 ,        1]  =  0.566659
  [6171412 ,        1]  =  -0.208408
  [6767957 ,        1]  =  -0.313653
  ⋮
  [5018886 ,  9999999]  =  -0.569702
  [5732979 ,  9999999]  =  0.338108
  [7234948 ,  9999999]  =  0.0112105
  [7356030 ,  9999999]  =  -0.2578
  [8047043 ,  9999999]  =  -0.204781
  [8304474 ,  9999999]  =  -0.666206
  [8677569 ,  9999999]  =  -1.24693
  [1103670 , 10000000]  =  -2.20066
  [1845980 , 10000000]  =  -1.25441
  [3347778 , 10000000]  =  -1.05728
  [3377766 , 10000000]  =  -1.67352
  [8828535 , 10000000]  =  -0.895992
In [28]:
##Note that getting a random column in S is fast
@time S[:,rand(1:size(T,2))]

  0.000011 seconds (8 allocations: 512 bytes)
Out[28]:
10000000-element SparseVector{Float64,Int64} with 7 stored entries:
  [591301  ]  =  1.42841
  [1432026 ]  =  -1.46346
  [1436064 ]  =  0.788249
  [5157199 ]  =  -1.11265
  [8375008 ]  =  0.995188
  [8421036 ]  =  -1.19796
  [9936245 ]  =  -0.205836
In [29]:
## But getting a random row in S is slow
@time S[rand(1:size(T,1)),:]

  0.076382 seconds (14 allocations: 928 bytes)
Out[29]:
10000000-element SparseVector{Float64,Int64} with 9 stored entries:
  [1277159 ]  =  0.95178
  [4346296 ]  =  0.699655
  [5358137 ]  =  1.36552
  [5359960 ]  =  -0.964119
  [5570648 ]  =  -0.504028
  [7128779 ]  =  1.5484
  [8033799 ]  =  -0.311346
  [8112310 ]  =  -0.168015
  [8848492 ]  =  -0.242558