# Gamma Process Non-negative Matrix Factorization (GaP-NMF)¶

This notebook demonstrates how GaP-NMF works on real data. The code used in the demo is a Julia translation of MATALB code for GaP-NMF as in:

Bayesian Nonparametric Matrix Factorization for Recorded Music by Matthew D. Hoffman et al. in ICML 2010

The original MATLAB code: http://www.cs.princeton.edu/~mdhoffma/code/gapnmfmatlab.tar

In [2]:
using PyPlot

INFO: Loading help data...

In [3]:
# Bayesian NMF package
using BNMF

In [4]:
# Read test file
using WAV
x = x[:] # monoral
println(length(x)/fs, " sec.")

3.79375 sec.

In [5]:
# Short-Time Fourier Transform
include("stft.jl")

Out[5]:
istft (generic function with 1 method)
In [6]:
# Compute spectrogram
framelen = 2048
X = stft(x, framelen=framelen, hopsize=256, window=hanning(framelen))
X = abs(X[:,1:framelen/2])

figure(figsize=(16, 6), dpi=80, facecolor="w", edgecolor="k")
imshow(log(X)', origin="lower", aspect="auto")
colorbar()

Out[6]:
PyObject <matplotlib.colorbar.Colorbar instance at 0x7f942299ab00>
In [7]:
p = GaPNMF(X, K=100, a=0.1, b=0.1, alpha=1.0, smoothness=100)
typeof(p)

Out[7]:
GaPNMF (constructor with 3 methods)
In [8]:
tic()
fit!(p, epochs=100, verbose=true)
toc()

#1 bound: -258355.66569987283,
improvement: NaN,
activek: 100
#2 bound: -166964.96002387482,
improvement: 0.35373989352401103,
activek: 100
#3 bound: -129991.8965370445,
improvement: 0.22144205276091128,
activek: 100
#4 bound: -114130.39875258185,
improvement: 0.12201912739954912,
activek: 100
#5 bound: -106144.52246008086,
improvement: 0.0699715096046691,
activek: 100
#6 bound: -100518.69946800527,
improvement: 0.05300153848439386,
activek: 100
#7 bound: -94373.87015499652,
improvement: 0.0611312058903491,
activek: 100
#8 bound: -85407.36706304806,
improvement: 0.09501044173797438,
activek: 100
#9 bound: -71379.20170161093,
improvement: 0.16425006230529837,
activek: 100
#10 bound: -51067.294076852355,
improvement: 0.284563390182888,
activek: 99
#11 bound: -24669.156508926593,
improvement: 0.5169284577365425,
activek: 94
#12 bound: 4102.394646412344,
improvement: 1.166296510581052,
activek: 86
#13 bound: 37172.20978047181,
improvement: 8.061100402171185,
activek: 80
#14 bound: 71905.2228691249,
improvement: 0.9343811759853963,
activek: 70
#15 bound: 104135.88179713141,
improvement: 0.4482380784309606,
activek: 59
#16 bound: 133533.96505598474,
improvement: 0.28230503023083,
activek: 46
#17 bound: 157475.44680815807,
improvement: 0.17929132668333303,
activek: 39
#18 bound: 176234.83052766905,
improvement: 0.11912576912618193,
activek: 33
#19 bound: 191371.97749562288,
improvement: 0.08589191434310302,
activek: 31
#20 bound: 203320.95952385763,
improvement: 0.0624385146906268,
activek: 28
#21 bound: 212527.40919362585,
improvement: 0.045280376855038086,
activek: 27
#22 bound: 221271.98254512594,
improvement: 0.041145626273236285,
activek: 24
#23 bound: 227047.40465476713,
improvement: 0.026101009459990537,
activek: 22
#24 bound: 231900.9935011162,
improvement: 0.021376984483610847,
activek: 22
#25 bound: 236935.8543801122,
improvement: 0.02171125187081945,
activek: 22
#26 bound: 242323.12157731538,
improvement: 0.022737239204668788,
activek: 19
#27 bound: 245692.53493397307,
improvement: 0.013904630044073822,
activek: 19
#28 bound: 249330.1082934242,
improvement: 0.014805388207780368,
activek: 19
#29 bound: 253188.99310633534,
improvement: 0.015477010936721,
activek: 19
#30 bound: 257098.1937642234,
improvement: 0.015439852301344915,
activek: 18
#31 bound: 260811.65862233413,
improvement: 0.014443760976073729,
activek: 18
#32 bound: 264488.21050513384,
improvement: 0.014096577975923657,
activek: 16
#33 bound: 267381.2119991927,
improvement: 0.010938111337868948,
activek: 15
#34 bound: 269830.27669571573,
improvement: 0.009159449454999127,
activek: 15
#35 bound: 272294.75787379453,
improvement: 0.009133449397370501,
activek: 15
#36 bound: 274581.35622993053,
improvement: 0.008397511483477829,
activek: 15
#37 bound: 276529.68533708015,
improvement: 0.00709563509300361,
activek: 14
#38 bound: 278124.43559486896,
improvement: 0.005767012882703217,
activek: 14
#39 bound: 279586.42900052277,
improvement: 0.005256616170840243,
activek: 14
#40 bound: 280897.39361979894,
improvement: 0.004688942249316833,
activek: 13
#41 bound: 281907.5343082436,
improvement: 0.0035961198337493686,
activek: 13
#42 bound: 282847.9552359242,
improvement: 0.003335919807848574,
activek: 13
#43 bound: 283730.3996598744,
improvement: 0.0031198543514805036,
activek: 13
#44 bound: 284563.1679619092,
improvement: 0.002935069005764316,
activek: 13
#45 bound: 285351.8363099131,
improvement: 0.0027715053696248074,
activek: 13
#46 bound: 286101.56184193055,
improvement: 0.002627372375495129,
activek: 13
#47 bound: 286818.95426385384,
improvement: 0.0025074746789381106,
activek: 13
#48 bound: 287506.9574985963,
improvement: 0.002398736989011982,
activek: 13
#49 bound: 288154.78159329866,
improvement: 0.002253246670406354,
activek: 13
#50 bound: 288743.0419019079,
improvement: 0.0020414733545511578,
activek: 12
#51 bound: 289252.8625989089,
improvement: 0.0017656553510099952,
activek: 12
#52 bound: 289725.232744448,
improvement: 0.0016330699073983674,
activek: 12
#53 bound: 290162.90757765586,
improvement: 0.0015106548679310501,
activek: 12
#54 bound: 290569.38496781344,
improvement: 0.0014008592399040452,
activek: 12
#55 bound: 290948.291487965,
improvement: 0.0013040139111472632,
activek: 12
#56 bound: 291303.1780167969,
improvement: 0.0012197580780314925,
activek: 12
#57 bound: 291636.2081829912,
improvement: 0.0011432424749417106,
activek: 12
#58 bound: 291948.6619722085,
improvement: 0.001071382017905087,
activek: 12
#59 bound: 292242.2464596381,
improvement: 0.0010056031270922305,
activek: 12
#60 bound: 292518.53416279465,
improvement: 0.0009454064444946287,
activek: 12
#61 bound: 292778.7895722327,
improvement: 0.0008897057076499888,
activek: 12
#62 bound: 293025.04199046985,
improvement: 0.0008410869468957971,
activek: 12
#63 bound: 293259.715969421,
improvement: 0.0008008666336401717,
activek: 12
#64 bound: 293484.2688095808,
improvement: 0.0007657132157326826,
activek: 12
#65 bound: 293699.2632813068,
improvement: 0.000732558758934513,
activek: 12
#66 bound: 293905.0015859551,
improvement: 0.0007005067099921809,
activek: 12
#67 bound: 294101.767618701,
improvement: 0.0006694885479461332,
activek: 12
#68 bound: 294290.2079882837,
improvement: 0.0006407318497556445,
activek: 12
#69 bound: 294471.29677259346,
improvement: 0.0006153408417753195,
activek: 12
#70 bound: 294645.7946388076,
improvement: 0.0005925802213208346,
activek: 12
#71 bound: 294813.6967608542,
improvement: 0.0005698439451764352,
activek: 12
#72 bound: 294974.629941982,
improvement: 0.0005458809509054502,
activek: 12
#73 bound: 295128.54492309236,
improvement: 0.0005217905727711023,
activek: 12
#74 bound: 295275.94797967165,
improvement: 0.000499453743512694,
activek: 12
#75 bound: 295417.893278301,
improvement: 0.00048072082944974404,
activek: 12
#76 bound: 295555.37588700384,
improvement: 0.000465383484991994,
activek: 12
#77 bound: 295688.6802731279,
improvement: 0.00045103015204513733,
activek: 12
#78 bound: 295817.46751446003,
improvement: 0.0004355501239112557,
activek: 12
#79 bound: 295941.1726987199,
improvement: 0.0004181807967570829,
activek: 12
#80 bound: 296059.4741081589,
improvement: 0.00039974636972681316,
activek: 12
#81 bound: 296172.4748717353,
improvement: 0.000381682646423729,
activek: 12
#82 bound: 296280.6788987094,
improvement: 0.000365341266169173,
activek: 12
#83 bound: 296384.9379170988,
improvement: 0.00035189273487872333,
activek: 12
#84 bound: 296485.85837681143,
improvement: 0.0003405046842861189,
activek: 12
#85 bound: 296583.4631691834,
improvement: 0.00032920555775016393,
activek: 12
#86 bound: 296677.33816635475,
improvement: 0.00031652134670027217,
activek: 12
#87 bound: 296766.91172191023,
improvement: 0.00030192247277498294,
activek: 12
#88 bound: 296852.27510867815,
improvement: 0.0002876445567082202,
activek: 12
#89 bound: 296934.38152826263,
improvement: 0.0002765901644325338,
activek: 12
#90 bound: 297014.2363602156,
improvement: 0.0002689309050099011,
activek: 12
#91 bound: 297092.3936934883,
improvement: 0.00026314339080350937,
activek: 12
#92 bound: 297169.241321607,
improvement: 0.0002586657543241824,
activek: 12
#93 bound: 297245.1230941944,
improvement: 0.00025534867690174215,
activek: 12
#94 bound: 297320.00587926933,
improvement: 0.00025192267006919736,
activek: 12
#95 bound: 297393.48991517304,
improvement: 0.0002471546967934055,
activek: 12
#96 bound: 297465.2228586288,
improvement: 0.00024120549335564673,
activek: 12
#97 bound: 297535.04699156235,
improvement: 0.0002347304073483458,
activek: 12
#98 bound: 297602.8387909525,
improvement: 0.000227844753334511,
activek: 12
#99 bound: 297668.60083128273,
improvement: 0.00022097249003866942,
activek: 12
#100 bound: 297732.4308089611,
improvement: 0.00021443302215992348,
activek: 12
iteration finished
elapsed time: 65.55415797 seconds

Out[8]:
65.55415797
In [9]:
# Get active K-components
good = goodk(p)

goodspec = p.Eh[good,:]

Out[9]:
12x1024 Array{Float64,2}:
0.00277722    0.00726524   0.0211036  â€¦  0.000520452  0.000491621
0.00280999    0.00493925   0.0108074     0.00037936   0.000340876
0.00271666    0.00594996   0.0126983     0.000312902  0.0166263
0.00580218    0.0131258    0.0379031     0.172967     0.125151
11.8068       34.5922      84.7486        0.000627676  0.000588681
0.00194749    0.0042438    3.31672    â€¦  0.249485     0.274738
3.31153       3.91104      4.46299       0.0971082    0.0279081
0.0036766     0.00851538   0.0226062     0.000484815  0.000428812
0.00377926    1.49767      4.74214       0.000353943  0.000278075
0.000794334   0.456527     1.37461       0.867147     0.964443
0.00158975    0.00809176   0.0213092  â€¦  0.446736     0.363154
0.464875      0.698002     0.84818       0.0435258    0.0458247  
In [10]:
# Visualize K-active spectrum
figure(figsize=(16, 6), dpi=80, facecolor="w", edgecolor="k")

for k=1:size(goodspec,1)
plot(goodspec[k,:][:], label=string("#", k))
end
legend()

Out[10]:
PyObject <matplotlib.legend.Legend object at 0x7f94225e8510>
In [11]:
# Reconstruction by goodk
Y = xbar(p, good)
figure(figsize=(16, 6), dpi=80, facecolor="w", edgecolor="k")
imshow(log(Y)', origin="lower", aspect="auto", interpolation="nearest")
colorbar()