Fun Q: QWAFAxNEW
Goals
- Introduce you to Q and KDB+
- Amaze you with 12 machine-learning algorithms
- Impress you with efficient code
- Dazzle you with ASCII plots
- Motivate you to learn more about ML and Q
- Cheer you up with a fun quote
Q Fundamentals
What is Q?
- The language behind the world’s fastest time series database
- Language Paradigms:
- Interpreted
- Dynamically Typed
- Vectorized
- Functional
- Event-driven
- Think Python but with native support for:
- SQL
- Parallelization
- Large Tables (> 1 Billion Rows)
- Distributed Computing
Why do I like Q?
- Personally manage extremely large datasets with a single machine
- Databases are directories
- Partitions are sub-directories
- Tables are sub-directories within partitions
- Columns are files within table directories
- Efficiently joins disparate datasets—making complex analysis easy
- Vectorized computations reduces the interpreter overhead
- Vectors, matrices, dictionaries and tables are all lists
- Notation works intuitively across types
- Code-injection allows iterative development
- Q-SQL permits all q operators and even user-defined functions
- Null and infinite values are available for integer and temporal types
- Functional paradigm promotes parallelism by separating functions from data
Fun Q Overview
Where does Fun Q Fit in?
- The original Coursera machine-learning course was taught in Matlab/Octave by Andrew Ng
- Why can’t q be used for the same problems?
- It turns out that many machine-learning algorithms have beautifully vectorized solutions in q
- 5 years of continuous refactoring revealed common functional paradigms in machine learning
- Small 1–5 line functions are the perfect tool for teaching
- 12 distinct machine-learning algorithm families are taught using nothing but the 32-bit ‘personal’ edition or 64-bit ‘on-demand’ edition of the q language
- Limiting all graphics to the terminal adds a bit of fun to the journey
Algorithms
- K-Nearest Neighbors (KNN)
- K-Means/Medians/Mediods Clustering
- Hierarchical Agglomerative Clustering (HAC)
- Expectation Maximization (EM)
- Naive Bayes
- Decision Tree (ID3,C4.5,CART)
- Discrete Adaptive Boosting (AdaBoost)
- Boosted Aggregating (BAG) and Random Forest
- Linear Regression
- Logistic Regression and One-vs.-All
- Neural Network Classification/Regression
- Content-Based/Collaborative Filtering (Recommender Systems)
- Google PageRank
K-Nearest Neighbors (KNN)
K-Nearest Neighbors (KNN) Algorithm
- Predictions are made by aggregating the K nearest neighbors
-
‘Nearness’ requires the definition of distance metrics
mnorm:sum abs:: / Manhattan (taxicab) norm enorm2:{x wsum x} / Euclidean norm squared enorm:sqrt enorm2:: / Euclidean norm hdist:sum (<>):: / Hamming distance mdist:mnorm (-):: / Manhattan (taxicab) distance edist2:enorm2 (-):: / Euclidean distance squared edist:enorm (-):: / Euclidean distance -
Use data type to differentiate between classification and regression
/ weighted average or mode isord:{type[x] in 0 8 9h} / is ordered aom:{$[isord x;avg;mode]x} / average or mode waom:{[w;x]$[isord x;nwavg;wmode][w;x]} / weighted average or mode
/ find (k) smallest values from (d)istance vector (or matrix) and use
/ (w)eighting (f)unction to return the best estimate of y
knn:{[wf;k;y;d]
if[not type d;:.z.s[wf;k;y] peach d]; / recurse for matrix d
if[any n:null d;d@:i:where not n; y@:i]; / filter null distances
p:(waom . (wf d::;y)@\:#[;iasc d]::) peach k&count d; / make predictions
p}
K-Nearest Neighbors (KNN) Example
-
Load Fun Q library and famous Iris data set
q)\l funq.q q)\l iris.q [down]loading iris data set q)iris.t species slength swidth plength pwidth ----------------------------------------- Iris-setosa 5.1 3.5 1.4 0.2 Iris-setosa 4.9 3 1.4 0.2 Iris-setosa 4.7 3.2 1.3 0.2 Iris-setosa 4.6 3.1 1.5 0.2 .. -
Partition between train/test (and even validation)
q)d:.ut.part[`test`train!4 1;0N?] iris.t q)y:first get first `Y`X set' 0 1 cut value flip d.train q)yt:first get first `Yt`Xt set' 0 1 cut value flip d.test -
Train and test
q)k:5 q)avg yt=.ml.knn[0n<;k;y] .ml.f2nd[.ml.edist X] Xt 0.975 q)avg yt=.ml.knn[sqrt 1f%;k;y] .ml.pedist2[X;Xt] 0.9583333
K-Means Clustering
K-Means Clustering Algorithm
- Find K centroids which minimize the intra-cluster distances
- Lloyd’s algorithm iteratively alternates between:
- Assigning each data point to its nearest centroid
- Updating centroid definitions by summarizing each cluster
- Fun Q algorithm implements a single round of Lloyd’s algorithm
- Iteration is controlled by the user using q adverbs
- Fixed number of iterations
- Run until convergence
- Run until custom test returns false
/ using the (d)istance (f)unction, group matri(X) based on the closest
/ (C)entroid and return the cluster indices
cgroup:{[df;X;C]value group imin f2nd[df X] C}
/ Stuart Lloyd's algorithm. uses (d)istance (f)unction to assign the
/ matri(X) to the nearest (C)entroid and then uses the (c)entroid (f)unction
/ to update the centroid location.
lloyd:{[df;cf;X;C]cf X@\: cgroup[df;X;C]}
kmeans:lloyd[edist2;avg''] / k-means
K-Means Clustering Example
-
Plot actual clusters of petal length/width
q).ut.plt (t.plength;t.pwidth;{distinct[x]?x} t.species) 3| " " | " @@@@ " | " @@@@ @" 2| " @@@@@@@" | " ++x# @ " | " +++++ @ " 1| " ++++ " | " . " | "... " 0| "... " -
Plot computed centroids within each cluster
q)C:.ml.kmeans[X] over last 3 .ml.kmeanspp[X]// 2#() q).ut.plt .ml.append[0N;X 2 3],'.ml.append[1] C 2 3 3| " " | " .... " | " .... ." 2| " ...@..." | " .... . " | " ..@.. . " 1| " .... " | " . " | "... " 0| ".@. "
Hierarchical Agglomerative Clustering (HAC)
Hierarchical Agglomerative Clustering (HAC) Algorithm
- Iteratively merges closest clusters — one at a time
- Only the distances from previous round are needed
-
Implements a single round of the Lance-Williams algorithm
lancewilliams:{[lf;D;a;L] n:count D; d:D@'di:imin peach D; / find closest distances if[null d@:i:imin d;:(D;a;L)]; j:di i; / find closest clusters c:$[9h=type lf;lf;lf(freq a)@/:(i;j;til n)]; / determine coefficients nd:sum c*nd,(d;abs(-/)nd:D (i;j)); / calc new distances D[;i]:D[i]:nd; / update distances D[;j]:D[j]:n#0n; / erase j a[where j=a]:i; / all elements in cluster j are now in i L:L,'(j;i); / append linkage stats (D;a;L)} -
The
.ml.linkfunction iterates and returns the linkage statisticslink:{[lf;D] D:@'[D;a:til count D;:;0n]; / define cluster assignments and ignore loops if[-11h=type lf;lf:get lf]; / dereference lf L:last .[lancewilliams[lf]] over (D;a;2#()); / obtain linkage stats L}
Lance-Williams Formula and Coefficients
\[ d_{(ij)k} = \alpha_i d_{ik} + \alpha_j d_{jk} + \beta d_{ij} + \gamma \lvert d_{ik} - d_{jk} \rvert \]
| Method | \(\alpha_i\) | \(\alpha_j\) | \(\beta\) | \(\gamma\) |
|---|---|---|---|---|
| Single | .5 | .5 | 0 | -.5 |
| Complete | .5 | .5 | 0 | .5 |
| Average | \(\frac{n_i}{n_i+n_j}\) | \(\frac{n_j}{n_i+n_j}\) | 0 | 0 |
| Weighted | .5 | .5 | 0 | 0 |
| Centroid | \(\frac{n_i}{n_i+n_j}\) | \(\frac{n_j}{n_i+n_j}\) | \(-\frac{n_i n_j}{(n_i+n_j)^2}\) | 0 |
| Median | .5 | .5 | -.25 | 0 |
| Ward | \(\frac{n_i + n_k}{n_i+n_j+n_k}\) | \(\frac{n_j + n_k}{n_i+n_j+n_k}\) | \(-\frac{n_k}{n_i+n_j+n_k}\) | 0 |
Hierarchical Agglomerative Clustering (HAC) Example
-
Build linkage statistics
q)X:iris.X q)L:.ml.link[`.ml.lw.ward] sqrt .ml.pedist2[X;X] -
Visually inspect clusters of different sizes
q)(,'/) (.ut.plt X[2 3],enlist .ut.ugrp@) each .ml.clust[L] 1+til 3 3| " " | " ++++ @@@@ +@@@ " | " ++++ + @@@@ @ @@@@ @" 2| " +++++++ @@@@@@@ +#@@@@@" | " ++++ + @@@@ @ ++++ @ " | " +++++ + @@@@@ @ +++++ @ " 1| " ++++ @@@@ ++++ " | " + . . " | "+++ ... ... " 0| "+++ ... ... "
Expectation Maximization (EM)
Expectation Maximization (EM) Algorithm
- Iteratively alternates between using the likelihood and MLE to:
- Assign a probability that points belongs to a distribution
- Compute new distribution parameters
- Using Gaussian distribution, can be considered a soft K-Means
- Introduces probability metrics to the toolbox
-
Implements a single round of Expectation Maximization
/ using (l)ikelihood (f)unction, (w)eighted (m)aximum likelihood estimator / (f)unction with prior probabilities (p)hi and distribution parameters / (THETA), optionally (f)fit (p)hi and perform expectation maximization em:{[fp;lf;wmf;X;phi;THETA] W:prb likelihood[0b;lf;X;phi;THETA]; / weights (responsibilities) if[fp;phi:avg each W]; / new phi estimates THETA:wmf[;X] peach W; / new THETA estimates (phi;THETA)}
Expectation Maximization (EM) Example
-
Initialized multi-variate Gaussian distribution parameters
q)`X`y set' iris`X`y; q)k:count distinct y q)phi:k#1f%k q)mu:X@\:/:neg[k]?count y q)SIGMA:k#enlist X cov\:/: X q)lf:.ml.gaussmvl q)mf:.ml.wgaussmvmle -
Fit Iris data to the distribution and compare clusters
q)pT:(.ml.em[1b;lf;mf;X]//) (phi;flip (mu;SIGMA)) q)p:.ml.imax .ml.likelihood[1b;.ml.gaussmvll;X] . pT q)m:.ml.mode each y group p q)avg y=m p 0.8866667
Confusion Matrix Example
-
A confusion matrix compares expected vs actual values
/ given true labels y and predicted labels p, return a confusion matrix cm:{[y;p] n:count u:asc distinct y,p; m:./[(n;n)#0;flip (u?p;u?y);1+]; t:([]y:u)!flip (`$string u)!m; t} -
Q tables clearly render the confusion matrix
q).ut.totals[`TOTAL] .ml.cm[y;m p] y | Iris-setosa Iris-versicolor Iris-virginica TOTAL ---------------| ------------------------------------------------ Iris-setosa | 50 0 0 50 Iris-versicolor| 0 49 1 50 Iris-virginica | 0 16 34 50 | 50 65 35 150
Naive Bayes
Naive Bayes Algorithm
- ‘Naively’ assumes conditional independence between features
- Does not require iterating
-
Uses weighted maximum likelihood estimator functions along with class labels to build a classifier
/ fit parameters given (w)eighted (m)aximization (f)unction returns a / dictionary with prior and conditional likelihoods fnb:{[wmf;w;y;X] if[(::)~w;w:count[y]#1f]; / handle unassigned weight pT:(odds g; (wmf . (w;X@\:) @\:) peach g:group y); pT} -
Classifies by using [log]likelihood function and picking class with highest probability density
/ using a [log](l)ikelihood (f)unction and prior probabilities (p)hi and / distribution parameters (T)HETA, perform naive Bayes classification pnb:{[l;lf;pT;X] d:{(x . z) y}[lf]'[X] peach pT[1]; / compute probability densities c:imax $[l;log[pT 0]+sum flip d;pT[0]*prd flip d]; c}
Naive Bayes Example
-
Use MLE to associate each class with prior and conditional likelihoods (not probabilities)
q)flip pT:.ml.fnb[.ml.wgaussmle/:;::;iris.y;iris.X] Iris-setosa | 0.3333333 (5.006 0.121764;3.428 0.140816;1.462 0.029556;0.246 0.0108.. Iris-versicolor| 0.3333333 (5.936 0.261104;2.77 0.0965;4.26 0.2164;1.326 0.038324) .. Iris-virginica | 0.3333333 (6.588 0.396256;2.974 0.101924;5.552 0.298496;2.026 0.0739.. -
Classification picks the class with the maximum likelihood density
q)avg iris.y=.ml.pnb[0b;.ml.gaussl;pT] iris.X .96
Decision Tree
Decision Tree Algorithm
- Recursively splits observations on features, progressively creating branches with more uniform target values
- Parameters include ‘gain functions’, ‘impurity functions’ and options to control tree depth
dt:{[cgf;ogf;ipf;opt;w;t]
if[(::)~w;w:n#1f%n:count t]; / compute default weight vector
if[1=count d:flip t;:(w;first d)]; / no features to test
opt:(`maxd`minss`minsl`ming`maxff!(0N;2;1;0;::)),opt; / default options
if[0=opt`maxd;:(w;first d)]; / check if we've reached max depth
if[identical a:first d;:(w;a)]; / check if all values are equal
if[opt[`minss]>count a;:(w;a)]; / check if insufficient samples
d:((neg floor opt[`maxff] count d)?key d)#d:1 _d; / sub-select features
d:{.[x isord z;y] z}[(cgf;ogf);(ipf;w;a)] peach d; / compute gains
d:(where (any opt[`minsl]>count each last::) each d) _ d; / filter on minsl
if[0=count d;:(w;a)]; / check if all leaves have < minsl samples
if[opt[`ming]>=first b:d bf:imax d[;0];:(w;a)]; / check gain of best feature
c:count k:key g:last b; / grab subtrees, feature names and count
/ distribute nulls down each branch with reduced weight
if[c>ni:null[k]?1b;w:@[w;n:g nk:k ni;%;c-1];g:(nk _g),\:n];
if[(::)~b 1;t:(1#bf)_t]; / don't reuse exhausted features
b[2]:.z.s[cgf;ogf;ipf;@[opt;`maxd;-;1]]'[w g;t g]; / split sub-trees
bf,1_b}
Decision Tree Impurity Functions
- Impurity functions distill the impurity of branches into a single value
-
Classification impurity functions are used for categorical values
misc:{1f-avg x=mode x} / misclassification gini:{1f-enorm2 odds group x} / Gini entropy:{neg sum x*log x:odds group x} / entropy - The Gini impurity function is used in the CART decision tree
algorithm (arguably because it is more efficient than the entropy
function which uses the
logoperator) -
Regression impurity functions are used for ordered values
mse:{enorm2[x-avg x]%count x} / mean squared error mae:{avg abs x-avg x} / mean absolute error
Decision Tree Projections
- Many famous decision tree implementations are just
parameterizations of
.ml.dt - Scikit-learn’s CART algorithm trees treats all features as ordered – requiring users to first us one-hot encoding
- Ross Quinlan’s C4.5 algorithm differentiates between categorical and ordered features
- Regression trees use the same algorithm, but use a mean squared error impurity function and predict values with the weighted average of the branch values instead of the weighted mode
aid:dt[sig;oig;wmse] / automatic interaction detection
thaid:dt[sig;oig;wmisc] / theta automatic interaction detection
id3:dt[ig;ig;wentropy] / iterative dichotomizer 3
q45:dt[gr;ogr;wentropy] / like c4.5
ct:dt[oig;oig;wgini] / classification tree
rt:dt[oig;oig;wmse] / regression tree
Decision Tree Example (ID3)
-
To play or not to play
q)\l weather.q q)weather.t Play Outlook Temperature Humidity Wind ----------------------------------------- No Sunny Hot High Weak No Sunny Hot High Strong Yes Overcast Hot High Weak Yes Rain Mild High Weak Yes Rain Cool Normal Weak No Rain Cool Normal Strong Yes Overcast Cool Normal Strong No Sunny Mild High Weak Yes Sunny Cool Normal Weak .. -
Categorical features cannot be reused
q)-1 .ml.ptree[0] .ml.id3[();::] weather.t; Yes (n = 14, err = 35.7%) | Outlook :: Overcast: Yes (n = 4, err = 0%) | Outlook :: Rain: Yes (n = 5, err = 40%) | | Wind :: Strong: No (n = 2, err = 0%) | | Wind :: Weak: Yes (n = 3, err = 0%) | Outlook :: Sunny: No (n = 5, err = 40%) | | Humidity :: High: No (n = 3, err = 0%) | | Humidity :: Normal: Yes (n = 2, err = 0%)
Decision Tree Example (CART)
- Each split is binary
- Features can be reused
q)d:.ut.part[`train`test!3 1;iris.t.species] iris.t
q)-1 .ml.ptree[0] tr:.ml.ct[();::] d`train;
Iris-setosa (n = 111, err = 66.7%)
| plength >[;2.45] 0: Iris-setosa (n = 37, err = 0%)
| plength >[;2.45] 1: Iris-versicolor (n = 74, err = 50%)
| | pwidth >[;1.7] 0: Iris-versicolor (n = 39, err = 7.7%)
| | | plength >[;4.95] 0: Iris-versicolor (n = 35, err = 0%)
| | | plength >[;4.95] 1: Iris-virginica (n = 4, err = 25%)
| | | | pwidth >[;1.55] 0: Iris-virginica (n = 3, err = 0%)
| | | | pwidth >[;1.55] 1: Iris-versicolor (n = 1, err = 0%)
| | pwidth >[;1.7] 1: Iris-virginica (n = 35, err = 2.9%)
| | | plength >[;4.85] 0: Iris-virginica (n = 3, err = 33.3%)
| | | | swidth >[;3.1] 0: Iris-virginica (n = 2, err = 0%)
| | | | swidth >[;3.1] 1: Iris-versicolor (n = 1, err = 0%)
| | | plength >[;4.85] 1: Iris-virginica (n = 32, err = 0%)
q)avg d.test.species=.ml.pdt[tr] d`test
0.9487179
Discrete Adaptive Boosting (AdaBoost)
Discrete Adaptive Boosting Algorithm
- An ensemble of weak learners generates a strong learner
- Weak learners merely need to classify with an accuracy > 50%
- Decision stumps are commonly used as weak learners
- Each additional learner adjusts weights to fit the observations incorrectly classified during the previous round
adaboost:{[tf;cf;w;t]
if[(::)~w;w:n#1f%n:count t]; / initialize weights
m:tf[w] t; / train model
p:cf[m] t; / make predictions
e:sum w*not p=y:first flip t; / compute weighted error
a:.5*log (c:1f-e)%e; / compute alpha (minimize exponential loss)
/ w*:exp neg a*y*p; / increase/decrease weights
/ w%:sum w; / normalize weights
w%:2f*?[y=p;c;e]; / increase/decrease and normalize weights
(m;a;w)}
Wisconsin Diagnostic Breast Cancer Data Set
-
Includes mean, max and standard deviation of tumor features
q)\l wdbc.q [down]loading wisconsin-diagnostic-breast-cancer data set q)10?wdbc.t diagnosis radius_mean texture_mean perimeter_mean area_mean smoothness_mean compactne.. -------------------------------------------------------------------------------------.. B 11.89 18.35 77.32 432.2 0.09363 0.1154 .. M 15.13 29.81 96.71 719.5 0.0832 0.04605 .. B 12.85 21.37 82.63 514.5 0.07551 0.08316 .. M 22.27 19.67 152.8 1509 0.1326 0.2768 .. B 9.465 21.01 60.11 269.4 0.1044 0.07773 .. B 9.683 19.34 61.05 285.7 0.08491 0.0503 .. B 11.54 14.44 74.65 402.9 0.09984 0.112 .. B 12.7 12.17 80.88 495 0.08785 0.05794 .. M 13.61 24.69 87.76 572.6 0.09258 0.07862 .. M 18.05 16.15 120.2 1006 0.1065 0.2146 .. -
Malignant vs. Benign diagnosis feature can be used for binomial classification
Discrete Adaptive Boosting Example
-
Training accuracy improves with extra components
q)d:.ut.part[`train`test!3 1;0N?] update -1 1 "M"=diagnosis from 11#/:wdbc.t q)stump:.ml.ct[(1#`maxd)!1#1]; k:50 q)m:.ml.fab[k;stump;.ml.pdt] d.train q).ut.plt avg d.train.diagnosis = .ml.pab[1+til k;.ml.pdt;m] d.train 1 | " " 0.98| " +" | " +++++++" 0.96| " ++++++ + " | " ++++++ + " | "+++ " 0.94| " + " 0.92| " " | "+ " 0.9 | " " -
Testing accuracy does as well, but not as consistently
q).ut.plt avg d.test.diagnosis = .ml.pab[1+til k;.ml.pdt;m] d.test 0.98| " + ++" | " ++ +++++ " 0.96| "++++++++ + ++++" | " ++++ + + " 0.94| " ++ " | " " | " " 0.92| "+ " | " " 0.9 | " "
Boosted Aggregating (BAG) and Random Forest
Boosted Aggregating Algorithm
- Instead of complicated pruning techniques, fit many independent trees and average the predictions
- Boosting involves random sampling (with replacement) of the data
- The difference between Boosted Aggregating and a Random Forest comes down to how the tree is grown, not the ensemble algorithm
- To increase tree diversity, a Random Forest adds artificial limits on the number of observed features at each tree level
/ generate (n) decision trees by applying (f) to a resampled (with
/ replacement) (t)able
bag:{[n;f;t](f ?[;t]::) peach n#count t} / (b)ootstrap (ag)gregating
/ given an atom or list (k), and bootstrap aggregating (m)odel, make
/ prediction on (d)ictionary
pbag:{[k;m;d]
if[count[m]<max k;'`length];
if[98h=type d;:.z.s[k;m] peach d]; / iterate on a table
p:k {aom x#y}\: pdt[;d] peach m;
p}
Random Forest Example
- The decision tree algorithm provides the
maxffoption which controls the maximum number of features visible at each branch -
Random Forest classification (
maxff = sqrt)q)k:20 q)d:.ut.part[`train`test!3 1;0N?] wdbc.t q)m:.ml.bag[k;.ml.q45[(1#`maxff)!1#sqrt;::]] d`train q)avg d.test.diagnosis=.ml.pbag[k;m] d`test 0.9370629 -
Random Forest regression (
maxff = 1%3)q)\l winequality.q [down]loading wine-quality data set q)d:.ut.part[`train`test!3 1;0N?] winequality.red.t q)m:.ml.bag[k;.ml.rt[(1#`maxff)!1#%[;3];::]] d`train q).ml.rms d.test.quality-.ml.pbag[k;m] d`test 0.5682814
Linear Regression
Linear Regression Algorithm
- Q has a native multi-variate
lsqoperator -
The implementation is similar to this normal equations function
/ given target matrix Y and data matri(X), return the THETA matrix resulting / from minimizing sum of squared residuals normeq:{[Y;X]mm[mmt[Y;X]] minv mmt[X;X]} / normal equations ols - Inverting matrices is slow and numerically unstable
- So the
lsqoperator uses Cholesky decomposition to solve for the coefficients directly instead of using LU decomposition to first invert the normal equations matrix. -
Ridge regression is simply the addition of L2 regularization
/ given (l2) regularization parameter, target matrix Y and data matri(X), / return the THETA matrix resulting from performing ridge regression ridge:{[l2;Y;X]mm[mmt[Y;X]] minv mmt[X;X]+diag count[X]#l2}
Linear Regression Example
-
Use Box-Muller to generate correlated random variables
q)plt:.ut.plot[30;15;.ut.c10] q)X:(.ml.bm 10000?) each 1 1f q)rho:.8 / correlation q)X[1]:(rho;sqrt 1f-rho*rho)$X q)plt[sum] X 4 | " . . " | " . . .. ." | " ............" 2 | " ............. " | " ....---::-..... " | " ....-:++++:-..... " | " ....:=+#%%x+--.. " 0 | " ....-=x%@@%x:-... " | " ....:+#%%%+:-... " | " ....-+++x=:-... . " | " ....---::--.... " -2| " .............. " | ". .......... " | " . ...... " -4| " .. " -
Prepending a vector of ones allows us to fit the intercept
q)show THETA:(-1#X) lsq .ml.prepend[1f] 1#X -0.008971004 0.8033666
Logistic Regression
Logistic Regression Cost Function
-
Logistic regression applies the Sigmoid function to compress values between 0 and 1: \(S(x)=\frac{1}{1+e^{-x}}\)
sigmoid:1f%1f+exp neg:: / sigmoid - Cost function with L1 and L2 regularization \[ J(\theta) = \frac{1}{m} \left[ \sum_{i=1}^m \text{LogLoss}(y^{(i)},h_\theta(x^{(i)})) \right] + \frac{\lambda_1}{m} \sum_{j=1}^n\left|\theta_j\right| + \frac{\lambda_2}{2m} \sum_{j=1}^n\theta_j^2\]
-
Regularized logistic regression cost function
/ logistic regression cost logcost:{[rf;Y;X;THETA] J:(1f%m:count X 0)*revo[sum] logloss[Y] plog[X;THETA]; / cost if[count rf,:();THETA[;0]:0f; J+:sum rf[;m][;0][;THETA]]; / regularization J}
Logistic Regression Gradient Function
- Gradient function with L1 and L2 regularization \[G(\theta_j)=\frac{\partial J(\theta)}{\partial \theta_j}= \frac{1}{m}\left(h_\theta(x^{(i)})-y^{(i)}\right)x_j^{(i)} + \frac{\lambda_1}{m} \text{sign}(\theta_j) + \frac{\lambda_2}{m} \theta_j \]
-
Regularized logistic regression gradient function
/ logistic regression gradient loggrad:{[rf;Y;X;THETA] G:(1f%m:count X 0)*mmt[sigmoid[mm[THETA;X]]-Y] X:prepend[1f] X; / gradient if[count rf,:();THETA[;0]:0f; G+:sum rf[;m][;1][;THETA]]; / regularization G}
Logistic Regression Example
-
Regularization choices have been factored out
l1:{[l;m]((l%m)*revo[sum] abs::;(l%m)*signum::)} l2:{[l;m]((.5*l%m)*revo[sum] {x*x}::;(l%m)*)} -
Elastic net is simply a combination of the two
enet:{[a;lr](l1 a*lr;l2 a*1f-lr)} -
Regularized gradient descent prevents over-fitting
q)t:11#/:update "M"=diagnosis from wdbc.t q)d:.ut.part[`train`test!3 1;0N?] "f"$t q)y:first get first `Y`X set' 0 1 cut value flip d`train q)yt:first get first `Yt`Xt set' 0 1 cut value flip d`test q)THETA:10000 .ml.gd[.1;.ml.loggrad[.ml.l2[.1];Y;X]]/ enlist (1+count X)#0f q)avg yt="i"$p:first .ml.plog[Xt] THETA 0.8601399 q).ut.totals[`TOTAL] .ml.cm . "i"$(yt;p) y| 0 1 TOTAL -| ----------- 0| 83 7 90 1| 13 40 53 | 96 47 143
Neural Network Classification/Regression
Neural Network Algorithm
- Conceptually similar to applying multiple rounds of linear and/or logistic regression in sequence
- But our choice of activation function and cost function can vary
- As can the number (and size of) of hidden layers
- Gradient descent is once again the workhorse
- Parameter adjustments are computed by applying the chain rule to the final error and intermediate gradients
/ (r)egularization (f)unction, (n)etwork topology dimensions, hgolf:
/ (h)idden (g)radient (o)utput (l)oss functions
nncostgrad:{[rf;n;hgolf;Y;X;theta]
THETA:nncut[n] theta; / unroll theta
ZA:enlist[(X;X)],(X;X) {(z;x z:plin[y 1;z])}[hgolf`h]\ -1_THETA;
P:hgolf[`o] plin[last[ZA]1;last THETA]; / prediction
J:(1f%m:count X 0)*revo[sum] hgolf[`l][Y;P]; / cost
G:hgolf[`g]@'`z`a!/:1_ZA; / activation gradient
D:reverse{[D;THETA;G]G*1_mtm[THETA;D]}\[E:P-Y;reverse 1_THETA;reverse G];
G:(1f%m)*(D,enlist E) mmt' prepend[1f] each ZA[;1]; / full gradient
if[count rf,:();THETA[;;0]:0f;JG:rf[;m][;;THETA];J+:sum JG@'0;G+:sum JG@'1];
(J;2 raze/ G)}
MNIST Hand Written Digits
-
Each digit is represented by a 28 x 28 pixel image with values ranging between 0 and 255
q)\l mnist.q [down]loading handwritten-digit data set q)`X`Xt`y`yt set' mnist`X`Xt`y`yt; q)X:1000#'X;y:1000#y; q)X%:255f;Xt%:255f -
The plotting utility allows us to visualize the images
q)plt:value .ut.plot[28;14;.ut.c10;avg] .ut.hmap flip 28 cut q)-1 (,'/) plt each X@\:/: -3?count X 0; :- .=#%@@x. . :=:=%@#-.-::x%: -%@% :x@%+--#+. =@@@@@@@@@@@@@@: .=: =@%- -+%#=. .:: -%@@#=--:---:x@@@- :%%- .#@= x@%- -:x@@%. .xx: :%@@@+. #@@=--...-@@: =@@%@@@@@x .x@@@x. +@@@@@@%%%@@@= :@@= .x@@@#- #@#:- :%@: :@@+ .x@@@#- =%@- +@@+ x@@@%- -#x. x@@: +@@@%: -:. -:::-
Neural Network Example
-
100 rounds of gradient descent
q)Y:.ml.diag[(1+max y)#1f]@\:y q)n:0N!"j"$.ut.nseq[2;count X;count Y] 784 397 10 q)theta:2 raze/ THETA:.ml.glorotu'[1+-1_n;1_n]; q)rf:.ml.l2[1f]; q)hgolf:`h`g`o`l!`.ml.sigmoid`.ml.dsigmoid`.ml.sigmoid`.ml.logloss q)theta:first .fmincg.fmincg[100;.ml.nncostgrad[rf;n;hgolf;Y;X];theta]; Iteration 100 | cost: 0.5257291 q)avg yt=p:.ml.imax .ml.pnn[hgolf;Xt] .ml.nncut[n] theta 0.8722 -
The confusion matrix is intuitive
q).ut.totals[`TOTAL] .ml.cm[yt;"i"$p] y| 0 1 2 3 4 5 6 7 8 9 TOTAL -| -------------------------------------------------- 0| 942 0 6 1 0 10 13 5 3 0 980 1| 0 1096 4 1 0 3 2 1 27 1 1135 2| 11 13 899 13 9 3 19 16 40 9 1032 3| 3 3 37 785 3 107 4 16 42 10 1010 4| 1 4 5 1 842 1 20 1 10 97 982 5| 13 5 10 30 20 730 18 12 34 20 892 6| 12 4 21 0 9 30 877 1 4 0 958 7| 4 17 30 13 8 7 1 915 0 33 1028 8| 21 14 6 17 15 43 21 16 793 28 974 9| 13 5 5 17 43 13 1 60 9 843 1009 | 1020 1161 1023 878 949 947 976 1043 962 1041 10000
Content-Based/Collaborative Filtering (Recommender Systems)
Generating Personal Movie Ratings
-
Loading the MovieLens data set
q)\l mlens.q [down]loading latest movielens data set "unzip -n ml-latest-small.zip" -
Append personal ratings
q)r:1!select `mlens.movie$movieId,rating:0n from mlens.movie q)r,:([]movieId:173 208 260 435 1197 2005 1968i;rating:.5 .5 4 .5 4 4 4f) q)r,:([]movieId:2918 4006 53996 69526 87520 112370i;rating:5 5 4 4 5 5f) q)select movieId,rating,movieId.title from r where not null rating movieId rating title ---------------------------------------------------- 173 0.5 "Judge Dredd" 208 0.5 "Waterworld" 260 4 "Star Wars: Episode IV - A New Hope" 435 0.5 "Coneheads" 1197 4 "Princess Bride, The" 1968 4 "Breakfast Club, The" 2005 4 "Goonies, The" 2918 5 "Ferris Bueller's Day Off" 4006 5 "Transformers: The Movie" 53996 4 "Transformers" 69526 4 "Transformers: Revenge of the Fallen" 87520 5 "Transformers: Dark of the Moon" 112370 5 "Transformers: Age of Extinction"
Content-Based Filtering Fitting
- Assume that each movie can be described by a list of boolean flags indicating which genre the movie is in
-
Building a matrix of these flags allows us to use linear regression
q)Y:value[r]1#`rating q)X:"f"$flip genre in/: value[mlens.movie]`genres q)theta:first 0N!THETA:(1;1+count X)#0f ,0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0f q)rf:.ml.l2[.1] q)theta:first .fmincg.fmincg[20;.ml.lincostgrad[rf;Y;X];theta] Iteration 20 | cost: 8.972662e-05 -
Theta coefficients describe my genre preferences
q){(5#x),-5#x}desc genre!1_theta 1980 | 1.772478 2010 | 1.546712 1970 | 0.6938292 Comedy | 0.408703 Sci-Fi | 0.4022778 Children| -0.2086186 Romance | -0.2254342 Drama | -0.4333146 Fantasy | -0.8020893 1990 | -2.581026
Content-Based Filtering Prediction
-
Scores are generated by multiplying my preferences by the genres
q)r:update score:first .ml.plin[X;enlist theta] from r -
Highest recommendations actually score higher than 5
q)select[5;>score] movieId,score,movieId.title from r movieId score title --------------------------------------------------- 85261 5.680565 "Mars Needs Moms" 92681 5.521147 "Journey 2: The Mysterious Island" 81564 5.442022 "Megamind" 4121 5.441443 "Innerspace" 8633 5.441443 "Last Starfighter, The" -
Highest recommendations ‘check all the boxes’
q)select[5;>score] movieId.genres from r genres ------------------------------------------------------------------- `genre$`Action`Adventure`Animation`Children`Comedy`Sci-Fi`IMAX`2010 `genre$`Action`Adventure`Comedy`Sci-Fi`IMAX`2010 `genre$`Action`Animation`Children`Comedy`Sci-Fi`IMAX`2010 `genre$`Action`Adventure`Comedy`Sci-Fi`1980 `genre$`Action`Adventure`Comedy`Sci-Fi`1980
User-User Collaborative Filtering
-
Build a user/movie matrix and demean each user to remove biases
q)n:20 q)m:exec distinct movieId from rating where n<(count;i) fby movieId q)R:value exec (movieId!rating) m by userId from rating where movieId in m q)R,:r[([]movieId:m);`rating] q)U:R-au:avg each R -
Recommendations are generated using KNN to average the ratings of the k ‘closest’ users
q)k:30 q)p:last[au]+.ml.fknn[1f-;.ml.cordist\:;k;U;0f^U] 0f^last U q)`score xdesc update score:p,movieId.title from ([]movieId:m)#r movieId| rating score title -------| ------------------------------------------- 1967 | 5.117033 "Labyrinth" 475 | 4.761649 "In the Name of the Father" 1276 | 4.470915 "Cool Hand Luke" 2067 | 4.470915 "Doctor Zhivago" 2968 | 4.421184 "Time Bandits" 2918 | 5 4.403138 "Ferris Bueller's Day Off" 908 | 4.346668 "North by Northwest" 71535 | 4.279887 "Zombieland" 8874 | 4.279887 "Shaun of the Dead"
Item-Item Collaborative Filtering
- User ratings are sparse and subject to frequent changes
- Item ratings have more support and are more stable
-
Item (instead of user) biases are removed
q)I-:ai:avg each I:flip R -
Recommendations are generated using KNN to average our own rating of each movie’s k most similar movies
q)D:((0^I) .ml.cosdist\:) peach 0^I q)p:ai+.ml.knn[1f-;k;last each I] D q)`score xdesc update score:p,movieId.title from ([]movieId:m)#r movieId| rating score title -------| ----------------------------------------------------------------- 912 | 5.39 "Casablanca" 926 | 5.379167 "All About Eve" 908 | 5.334211 "North by Northwest" 356 | 5.314134 "Forrest Gump" 1262 | 5.277907 "Great Escape, The" 1207 | 5.27069 "To Kill a Mockingbird" 720 | 5.242593 "Wallace & Gromit: The Best of Aardman Animation" 899 | 5.224468 "Singin' in the Rain" 1247 | 5.213291 "Graduate, The" 1148 | 5.185714 "Wallace & Gromit: The Wrong Trousers"
Matrix Decomposition Collaborative Filtering Algorithm
- The null-riddled rating matrix cannot be decomposed with SVD
- Alternatives include Gradient Descent, Stochastic Gradient Descent and Alternating Least Squares (ALS).
- ALS elegantly transforms a non-convex problem into an iterative quadratic problem
/ ALS-WR (a)lternating (l)east (s)quares with (w)eighted (r)egularization
alswr:{[l2;Y;XTHETA]
X:flip f2nd[wridge[l2;;XTHETA 1]] Y; / hold THETA constant, solve for X
THETA:flip wridge[l2;;X] peach Y; / hold X constant, solve for THETA
(X;THETA)}
Matrix Decomposition Collaborative Filtering Example
-
Perform regularized ALS until cost is below supplied threshold
q)nf:10; q)n:(ni:count U 0;nu:count U) q)XTHETA:(X:-1+ni?/:nf#1f;THETA:-1+nu?/:nf#2f) q)XTHETA:first .ml.iter[1;.0001;.ml.cfcost[();U] .;.ml.alswr[.01;U]] XTHETA iter: 144 | cost: 8.349892 | pct: 9.952313e-05 -
Recommendations are generated by combining item and user exposure matrices
q)P:au+.ml.pcf . XTHETA q)show t:`score xdesc update score:last P,movieId.title from ([]movieId:m)#r movieId| rating score title -------| --------------------------------------------------------------------- 4232 | 9.907802 "Spy Kids" 102903 | 9.573078 "Now You See Me" 3869 | 8.985205 "Naked Gun 2 1/2: The Smell of Fear, The" 265 | 8.948548 "Like Water for Chocolate (Como agua para chocolate)" 84954 | 8.573156 "Adjustment Bureau, The" 116823 | 8.445161 "The Hunger Games: Mockingjay - Part 1" 53121 | 8.140992 "Shrek the Third" 6942 | 7.996043 "Love Actually"
Google PageRank
Google PageRank Algorithm
- The PageRank algorithm relies on page links (not content)
- Links (like citations) provide evidence of relative importance
-
The PageRank algorithm adds ‘random surfing’
/ given a (d)amping factor (1 - the probability of random surfing) and the / (A)djacency matrix, create the Markov Google matrix google:{[d;A] M:A%1f|s:sum each A; / convert to Markov matrix M+:(0f=s)%n:count M; / add links to dangling pages M:(d*M)+(1f-d)%n; / dampen M} -
The Google Matrix is a transition matrix with random surfing
q).ml.google[.85] (10110b;11001b;01110b;00111b;00000b) 0.3133333 0.03 0.3133333 0.3133333 0.03 0.3133333 0.3133333 0.03 0.03 0.3133333 0.03 0.3133333 0.3133333 0.3133333 0.03 0.03 0.03 0.3133333 0.3133333 0.3133333 0.2 0.2 0.2 0.2 0.2
Google PageRank Example
-
Generate page links
q)i:1 1 2 2 3 3 3 4 6 q)j:2 6 3 4 4 5 6 1 1 -
Enumerate links and create a sparse matrix
q)node:asc distinct raze l q)l:node?l q)show S:(1 2#1+max over l), .ml.prepend[1f] l 6 6 1 1 1 1 1 1 1 1 1f 0 0 1 1 2 2 2 3 5 1 5 2 3 3 4 5 0 0 -
Iteratively apply transition matrix until convergence
q)node[i]!r i:idesc r:.ml.pageranks[d;S] over r:n#1f%n:S[0;0] 1| 0.3210169 6| 0.200744 2| 0.170543 4| 0.1367926 3| 0.1065916 5| 0.0643118
Sparklines
Sparkline Algorithm
A sparkline is a small intense, simple, word-sized graphic with typographic resolution. – Edward Tufte
- Bucket all values into 8 bins
- Map bins to 8 Unicode characters
/ allocate x into n bins
nbin:{[n;x](n-1)&floor n*.5^x%max x-:min x}
/ generate unicode sparkline
spark:raze("c"$226 150,/:129+til 8)nbin[8]::
Sparkline Example
-
Load sample DJIA data set
q)\l dji.q [down]loading dji data set q)dji.t quarter stock date open high low close volume percent_change_price percent_change_vo.. -------------------------------------------------------------------------------------------------.. 1 AA 2011.01.07 15.82 16.72 15.78 16.42 239655616 3.79267 .. 1 AA 2011.01.14 16.71 16.71 15.64 15.97 242963398 -4.42849 1.380223 .. 1 AA 2011.01.21 16.19 16.38 15.6 15.79 138428495 -2.47066 -43.02496 .. 1 AA 2011.01.28 15.87 16.63 15.82 16.13 151379173 1.63831 9.3555 .. 1 AA 2011.02.04 16.18 17.39 16.18 17.14 154387761 5.93325 1.987452 .. .. -
Use q-sql to conveniently access
closeandstockcolumnsq)-1@'10#exec ((4$string stock 0),": ",.ut.spark close) by stock from dji.t; AA : ▅▄▃▄▇▇▇▅▅▄▄▆▇█▅▆▆▇▆▄▅▃▂▁▂ AXP : ▁▃▃▁▁▄▃▁▁▁▁▃▂▃▃▄▆▇▇██▇▅▆▆ BA : ▁▁▂▁▂▃▃▃▃▂▁▄▄▄▃▅███▇▆▅▃▄▂ BAC : ▇█▇▆▇██▇▇▇▆▅▅▆▄▄▃▄▃▂▂▂▁▁▁ CAT : ▁▁▁▂▃▄▅▄▄▃▅▆█▇▆▆█▇▅▅▅▃▂▂▃ CSCO: ▇█▇▇█▅▅▅▄▄▃▃▃▄▃▃▃▃▃▂▂▂▁▁▁ CVX : ▁▁▂▁▃▃▄▅▆▄▆▇██▇██▆▅▅▆▅▄▄▃ DD : ▂▂▁▂▄▆█▆▆▅▅▆▇▇▇██▆▅▅▄▂▂▂▄ DIS : ▃▃▃▂▅█████▅██▆▆▇██▆▆▆▃▂▁▁ GE : ▂▂▅▆▆██▇▆▆▃▅▆▆▅▅▆▅▅▄▄▂▁▂▁
Thank You
Never, ever underestimate the importance of having fun – Randy Pausch