Kernels
Distance and kernel matrices, row-wise argmin and argmax masks, kernel centring and a numerically stable softmax, for writing learners over.
This module provides the matrices that distance-based and kernel-based learners are written over: the Euclidean distance matrix between two sets of points, the linear, polynomial, RBF, Laplacian and sigmoid kernels, kernel centring, and a softmax and log-sum-exp that stay finite on large scores. It also answers a question a contraction cannot ask directly — which entry of each row is smallest or largest — by returning a mask of ones and zeros that marks it.
Conventions. Points are the rows of a matrix. x is n×p (n points) and c is K×p (K points or
centres, with the same p columns); every pairwise result is n×K, with row i holding point i against
every row of c. Per-row results are vectors n. Kernel parameters — the scale g, offset r and
degree d — are scalars.
A mask marks every position that attains the row’s minimum or maximum, so ties mark several positions.
That is deliberate: the data supports each of them equally, and choosing one would be arbitrary. Use
markedRow to read a value out of a mask; it averages over ties.
Broadcasting aligns on the trailing axis, so a vector n added to an n×K matrix would line up against
the K columns instead of the n rows and silently give a wrong answer. The definitions here lift every
per-row vector over the columns explicitly with lifted; do the same in your own formulas.
open Prelude.Kernels
Vectors of ones and shaped zeros
zeroed
zeroed(a) · a tensor → zero tensor shaped like a
Let \(A\) (the argument a) be a tensor of any shape. Then \(\operatorname{zeroed}(A)\) has the shape of \(A\) and entries
A tensor of zeros shaped like a, written so that no law in the Prelude can reduce it to an unshaped 0.
The obvious spelling a − a is reduced by Arithmetic’s cancels law in any document where that module
is open, which loses the shape; use zeroed wherever a shaped zero matters, such as the starting
coefficients of an iterative fit.
Intuition Zero everywhere, but in a form the rewrite laws cannot see through. a − a looks like zero to Arithmetic’s cancels and loses its shape; min(a^0, 0) is also zero everywhere, yet matches no law, so the shape survives into the generated code.
Use it when an iterative fit or loop needs a starting vector or matrix of zeros that will later be read by a contraction — initial coefficients, accumulators, starting weights. Where the zero is only ever added to a shaped value, LinearAlgebra’s zeros is fine; the bug zeroed prevents only appears in documents where Arithmetic is open.
shape-handlingoptimisation
def zeroed(a) ≔ min(ones(a), 0)
Building blocks
rowOnes(a)— Let \(A \in \mathbb{R}^{n \times p}\). Then \(\operatorname{rowOnes}(A)_i = \frac{1}{p} \sum_{j=1}^{p} A_{ij}^0 = 1\) for \(i = 1, \dots, n\), that is \(\operatorname{rowOnes}(A) = \mathbf{1}_n\). A vector of ones as long as the rows ofa(n for an n×p matrix).colOnes(a)— Let \(A \in \mathbb{R}^{n \times p}\). Then \(\operatorname{colOnes}(A)_j = \frac{1}{n} \sum_{i=1}^{n} A_{ij}^0 = 1\) for \(j = 1, \dots, p\), that is \(\operatorname{colOnes}(A) = \mathbf{1}_p\). A vector of ones as long as the columns ofa(p for an n×p matrix).
Row-wise extremes and masks
rowLeast
rowLeast(d) · d matrix n×K → vector n
Let \(D \in \mathbb{R}^{n \times K}\) (the argument d). Then
The smallest entry of each row.
Intuition For each point, its best (lowest) score across all the candidates in its row — the distance to its nearest centre, the cost of its cheapest option.
Use it when you need the value, not the position: the k-means objective (sum of each point’s squared distance to its nearest centre), a nearest-neighbour distance for outlier scoring, or a minimum cost per state. When you need which column attains it, use lowest.
clusteringoutlier-detectiondata-summary
def rowLeast(d) ≔ minimum(d, 1)
rowMost
rowMost(d) · d matrix n×K → vector n
Let \(D \in \mathbb{R}^{n \times K}\) (the argument d). Then
The largest entry of each row.
Intuition For each row, the best score on offer — the value of the best action in a state, the confidence of the most likely class, the strongest similarity.
Use it when you build a Bellman target (\(r + \gamma \max_a Q(s', a)\) with one row per state), report the top class probability, or shift scores for a stable softmax (as steadied does). When you need which column attains the maximum, use highest.
value-learningclassificationdecision-making
def rowMost(d) ≔ maximum(d, 1)
lifted
lifted(d, v) · d matrix n×K · v vector n → matrix n×K
Let \(D \in \mathbb{R}^{n \times K}\) (the argument d) and \(v \in \mathbb{R}^n\). Then
The per-row vector v repeated across the K columns of d, so that it can be compared with or combined
with d entry by entry.
Intuition Copy one number per row across all the columns, so a per-row value lines up with the matrix it belongs to. Broadcasting would instead align the vector with the last axis — the columns — and give a silently wrong answer whenever the sizes happen to fit.
Use it when you compare a matrix with a per-row quantity (entries equal to the row minimum, as lowest does), subtract a per-row mean or maximum, or divide by per-row totals (as softmaxRows does). Only the shape of d is used. Reach for it whenever a vector n meets a matrix n×K.
shape-handlingfeature-scaling
def lifted(d, v) ≔ einsum("i,j->ij", v, colOnes(d))
lowest
lowest(d) · d matrix n×K → mask n×K
Let \(D \in \mathbb{R}^{n \times K}\) (the argument d). Then
Hilbert returns this mask, not the index of the minimum.
A mask with 1 where an entry equals its row’s minimum and 0 elsewhere: a row-wise argmin. Applied to a distance matrix, it marks each point’s nearest centre, which is what k-means, nearest neighbours, decision trees and Viterbi decoding need. Ties mark every tied position.
Intuition “Which one wins?” answered as a one-hot table rather than an index. A graph of contractions cannot look up position \(j\), but it can multiply by a 0/1 mask, so the mask is the form of argmin the rest of a formula can use.
Use it when you assign points to their nearest centre (k-means, nearest-centroid classification), pick the cheapest option per row, or select matching values from another n×K matrix by multiplying. Beware of ties: a row can carry several ones, so counting or summing through the mask double-counts; markedRow averages them. For the minimum value alone, rowLeast is cheaper.
clusteringclassificationmasking
def lowest(d) ≔ d = lifted(d, rowLeast(d))
highest
highest(d) · d matrix n×K → mask n×K
Let \(D \in \mathbb{R}^{n \times K}\) (the argument d). Then
A mask with 1 where an entry equals its row’s maximum: a row-wise argmax. Ties mark every tied position.
Intuition The argmax as a one-hot table: 1 where a row reaches its maximum. Multiplying by it keeps only the winning entries, which is how “take the best” is written without an index.
Use it when you turn class scores or probabilities into predicted labels, choose the greedy action from a table of action values, or keep the best predecessor in max-product decoding. Ties mark several columns; break them yourself (for example with a tiny caller-supplied random offset) if you need exactly one. For the maximum value alone use rowMost.
classificationdecision-makingmasking
def highest(d) ≔ d = lifted(d, rowMost(d))
markedRow
markedRow(d, m) · d matrix n×K · m mask n×K → vector n
Let \(D \in \mathbb{R}^{n \times K}\) (the argument d) and \(M \in \{0, 1\}^{n \times K}\) (the argument m). Then
where \(s_i = \sum_{j=1}^{K} M_{ij}\) is the number of marks in row \(i\).
For each row, the value of d at the positions marked in m, averaged when several are marked. A row
with no mark gives 0.
Intuition Read one value per row through a mask: the mask says where to look, and the division by the number of marks turns several tied marks into their average instead of their sum.
Use it when the mask comes from one matrix and the values from another of the same shape — the distance to the centre chosen by lowest, the reward of the action chosen by highest, the label score of the winning class. A row with no mark reads 0, which is indistinguishable from a real zero, so make sure every row has at least one mark.
maskingclusteringclassification
def markedRow(d, m) ≔ einsum("ij,ij->i", d, m) / safe(einsum("ij->i", m))
Building blocks
spanned(d)— Let \(D \in \mathbb{R}^{n \times K}\). Then \(\operatorname{spanned}(D) \in \mathbb{R}^{n \times K \times K}\) with \(\operatorname{spanned}(D)_{ijk} = D_{ij}\). An n×K×K tensor holding d[i,j] at every k.facing(d)— Let \(D \in \mathbb{R}^{n \times K}\). Then \(\operatorname{facing}(D) \in \mathbb{R}^{n \times K \times K}\) with \(\operatorname{facing}(D)_{ijk} = D_{ik}\). An n×K×K tensor holding d[i,k] at every j.slowLowest(d)— Let \(D \in \mathbb{R}^{n \times K}\). Then \(\operatorname{slowLowest}(D)_{ij} = \big[\, \sum_{k=1}^{K} [\,D_{ij} \le D_{ik}\,] = K \,\big]\), which equals \(\operatorname{lowest}(D)_{ij}\). The same mask aslowest, computed exactly by comparing every entry with every other in its row; a reference form that costs one index more.slowHighest(d)— Let \(D \in \mathbb{R}^{n \times K}\). Then \(\operatorname{slowHighest}(D)_{ij} = \big[\, \sum_{k=1}^{K} [\,D_{ij} \ge D_{ik}\,] = K \,\big]\), which equals \(\operatorname{highest}(D)_{ij}\). The same mask ashighest, by the same exhaustive comparison.
The distance matrix
The squared Euclidean distance between every point of x and every point of c is computed from norms and
one cross product rather than from an n×K×p tensor of differences, which would be p times larger:
squared
squared(x, c) · x points n×p · c points K×p → matrix n×K
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) have rows \(x_i^{\top}\) and \(C \in \mathbb{R}^{K \times p}\) (the argument c) have rows \(c_k^{\top}\). Then
The squared Euclidean distances \(\lVert x_i - c_k \rVert^2\). Because the three terms nearly cancel for a point against itself, rounding can leave a tiny negative value.
Intuition The squared gap between every point and every centre, obtained from lengths and dot products instead of subtracting every pair. That keeps memory at n×K rather than n×K×p, at the price that nearly identical points can come out as tiny negative numbers.
Use it when you rank neighbours or assign clusters — the square root does not change the order, so skip it — or feed an RBF kernel or a k-means objective. Cost is \(O(nKp)\). When the actual distance matters or a negative would be harmful (a square root, a log), use distance, which floors at zero.
similarity-searchclusteringgeometry
def squared(x, c) ≔ heldNorms(x, c) + laidNorms(x, c) − 2 · crossTerm(x, c)
distance
distance(x, c) · x points n×p · c points K×p → matrix n×K
Let \(X \in \mathbb{R}^{n \times p}\) and \(C \in \mathbb{R}^{K \times p}\) be as in squared, and \(S = \operatorname{squared}(X, C)\). Then
which is \(\lVert x_i - c_k \rVert\) up to rounding.
The Euclidean distances \(\lVert x_i - c_k \rVert\). Squared distances are floored at zero before the
square root, so a point against itself gives exactly 0 rather than NaN.
Intuition The straight-line distance between every pair of points, in the units of the data. The floor at zero absorbs the rounding that makes a point’s distance to itself come out slightly negative, so the root never produces NaN.
Use it when you need real distances: k-nearest neighbours with distance weighting, medoids, silhouette scores, radius queries, or the Laplacian kernel. Features on different scales dominate unequally, so standardise columns first. To rank or compare, squared gives the same order without the square root.
similarity-searchclusteringoutlier-detection
def distance(x, c) ≔ sqrt(max(squared(x, c), 0))
Building blocks
selfNorms(x)— Let \(X \in \mathbb{R}^{n \times p}\) with rows \(x_i^{\top}\). Then \(\operatorname{selfNorms}(X)_i = \sum_{j=1}^{p} X_{ij}^2 = \lVert x_i \rVert^2\). The squared norm of each row ofx, a vector n.heldNorms(x, c)— Let \(X \in \mathbb{R}^{n \times p}\) and \(C \in \mathbb{R}^{K \times p}\). Then \(\operatorname{heldNorms}(X, C)_{ik} = \lVert x_i \rVert^2\) for \(k = 1, \dots, K\), that is \(\operatorname{selfNorms}(X) \, \mathbf{1}_K^{\top}\). The squared norms ofxspread into an n×K matrix.laidNorms(x, c)— Let \(X \in \mathbb{R}^{n \times p}\) and \(C \in \mathbb{R}^{K \times p}\) with rows \(c_k^{\top}\). Then \(\operatorname{laidNorms}(X, C)_{ik} = \lVert c_k \rVert^2\), that is \(\mathbf{1}_n \operatorname{selfNorms}(C)^{\top}\). The squared norms ofcspread into an n×K matrix.crossTerm(x, c)— Let \(X \in \mathbb{R}^{n \times p}\) and \(C \in \mathbb{R}^{K \times p}\). Then \(\operatorname{crossTerm}(X, C)_{ik} = \sum_{j=1}^{p} X_{ij} C_{kj} = x_i^{\top} c_k\), that is \(X C^{\top}\). The inner products of every row ofxwith every row ofc.
Kernels
Each kernel compares every row of x with every row of c and returns an n×K matrix. Passing the same
matrix twice gives the square Gram matrix a kernel method needs.
linear
linear(x, c) · x points n×p · c points K×p → matrix n×K
Let \(X \in \mathbb{R}^{n \times p}\) (the argument x) have rows \(x_i^{\top}\) and \(C \in \mathbb{R}^{K \times p}\) (the argument c) have rows \(c_k^{\top}\). Then
The linear kernel.
Intuition Similarity as the plain dot product: large when two points point the same way and are long. A kernel method run on it is just the ordinary linear method, written in terms of pairs of points instead of features.
Use it when you want a baseline before trying a nonlinear kernel, when features already number in the thousands (text, one-hot data) and a linear boundary is enough, or when rows are unit-normalised and the dot product is the cosine similarity. It is the cheapest kernel, \(O(nKp)\), and the building block of polynomial and sigmoidal.
similarity-searchclassificationregression
def linear(x, c) ≔ einsum("ij,kj->ik", x, c)
polynomial
polynomial(x, c, g, r, d) · x points n×p · c points K×p · g scale, scalar · r offset, scalar · d degree, scalar → matrix n×K
Let \(X \in \mathbb{R}^{n \times p}\) and \(C \in \mathbb{R}^{K \times p}\) be as in linear, and let \(g, r, d \in \mathbb{R}\). Then
The polynomial kernel.
Intuition A dot product raised to a power equals a dot product over all products of up to \(d\) features, so a linear method on this kernel learns feature interactions without ever building them. The offset \(r > 0\) keeps the lower-order terms; \(r = 0\) keeps only products of exactly degree \(d\).
Use it when interactions between features matter — degree 2 captures every pairwise interaction, as in many text and image tasks. Use a whole-number degree (2 or 3 is typical): a fractional power of a negative base is NaN. Values grow like the \(d\)-th power of the dot product, so scale features or shrink \(g\). For smooth, local similarity, rbf is usually the better default.
classificationregressionfeature-interactions
def polynomial(x, c, g, r, d) ≔ (g · linear(x, c) + r)^d
rbf
rbf(x, c, g) · x points n×p · c points K×p · g scale, scalar → matrix n×K
Let \(X \in \mathbb{R}^{n \times p}\) and \(C \in \mathbb{R}^{K \times p}\) be as in linear, let \(g \in \mathbb{R}\), and let \(S = \operatorname{squared}(X, C)\). Then
The Gaussian (radial basis function) kernel. The squared distance is used as computed, without the floor at zero that distance applies.
Intuition A bell-shaped similarity: 1 for identical points, falling smoothly toward 0 with distance. \(g\) sets the reach — the similarity halves at a distance of \(\sqrt{\ln 2 / g}\) — and \(g = 1/(2\sigma^2)\) turns it into an unnormalised Gaussian of width \(\sigma\).
Use it when you need a general-purpose nonlinear kernel for SVMs, kernel ridge, Gaussian processes or kernel PCA, or Gaussian bumps for density estimation and smoothing. Start \(g\) near \(1/(2 m^2)\) with \(m\) the median pairwise distance, then tune. Too large a \(g\) makes the matrix nearly the identity (memorising); too small makes every entry near 1. Standardise features first.
classificationregressiondensity-estimation
def rbf(x, c, g) ≔ exp(0 − g · squared(x, c))
laplacian
laplacian(x, c, g) · x points n×p · c points K×p · g scale, scalar → matrix n×K
Let \(X \in \mathbb{R}^{n \times p}\) and \(C \in \mathbb{R}^{K \times p}\) be as in linear, let \(g \in \mathbb{R}\), and let \(D = \operatorname{distance}(X, C)\). Then
The Laplacian kernel.
Intuition Like rbf, but decaying with the distance rather than its square: a sharp peak at zero and long tails. Nearby points are told apart more strongly and far points still keep some similarity, so the fitted function can have kinks rather than being perfectly smooth.
Use it when the target is rough or piecewise, or when an RBF’s similarities die off too fast for the spacing of your data. \(g\) is an inverse length scale: similarity falls to \(1/e\) at distance \(1/g\). It needs a square root per entry, so it costs a little more than rbf.
classificationregressionsimilarity-search
def laplacian(x, c, g) ≔ exp(0 − g · distance(x, c))
sigmoidal
sigmoidal(x, c, g, r) · x points n×p · c points K×p · g scale, scalar · r offset, scalar → matrix n×K
Let \(X \in \mathbb{R}^{n \times p}\) and \(C \in \mathbb{R}^{K \times p}\) be as in linear, and let \(g, r \in \mathbb{R}\). Then
The sigmoid (hyperbolic tangent) kernel. It is not positive
semi-definite for every choice of g and r.
Intuition A dot product squashed into \((-1, 1)\), like the activation of a single neural unit. An SVM using it behaves like a two-layer neural network whose hidden units sit at the training points.
Use it when you want to compare against a neural-network-style model or reproduce a method that specifies it. Because it is not positive semi-definite for every g and r, solvers that assume a valid kernel (Cholesky, conjugate gradient, kernel PCA) can fail or give meaningless results; a small positive \(g\) and negative \(r\) are the usual choice. For most problems rbf is the safer default.
classificationsimilarity-search
def sigmoidal(x, c, g, r) ≔ tanh(g · linear(x, c) + r)
centreKernel
centreKernel(k) · k kernel matrix n×n → matrix n×n
Let \(K \in \mathbb{R}^{n \times n}\) (the argument k). Then
or, in matrix form, \(\tilde{K} = K - \tfrac{1}{n} K \mathbf{1}_n \mathbf{1}_n^{\top} - \tfrac{1}{n} \mathbf{1}_n \mathbf{1}_n^{\top} K + \tfrac{1}{n^2} \mathbf{1}_n \mathbf{1}_n^{\top} K \mathbf{1}_n \mathbf{1}_n^{\top}\).
The centred kernel matrix, equivalent to centring the points in feature space without ever computing the feature map — the step kernel PCA rests on. Each term is a row mean, a column mean or the overall mean, spread back over the matrix.
Intuition PCA needs centred data, but a kernel’s features are never built, so their mean cannot be subtracted directly. Removing row, column and overall means from the kernel matrix gives exactly the inner products the centred features would have had.
Use it when you do kernel PCA (centre, then take eigenvectors), or compare kernel matrices with a measure that assumes centring, such as kernel alignment or HSIC-style dependence tests. It is meant for the square Gram matrix of one set of points; centring a train-versus-test kernel requires the training means, which this function does not use.
dimensionality-reductionhypothesis-testing
def centreKernel(k) ≔ k − kernelDown(k) − kernelAcross(k) + kernelMean(k)
Building blocks
kernelRowMeans(k)— Let \(K \in \mathbb{R}^{n \times m}\). Then \(\operatorname{kernelRowMeans}(K)_i = \frac{1}{m} \sum_{j=1}^{m} K_{ij}\) for \(i = 1, \dots, n\). The mean of each row ofk.kernelColMeans(k)— Let \(K \in \mathbb{R}^{n \times m}\). Then \(\operatorname{kernelColMeans}(K)_j = \frac{1}{n} \sum_{i=1}^{n} K_{ij}\) for \(j = 1, \dots, m\). The mean of each column ofk.kernelMean(k)— Let \(K \in \mathbb{R}^{n \times m}\). Then \(\operatorname{kernelMean}(K) = \frac{1}{nm} \sum_{i=1}^{n} \sum_{j=1}^{m} K_{ij}\). The mean of every entry ofk.kernelDown(k)— Let \(K \in \mathbb{R}^{n \times m}\) and \(\rho = \operatorname{kernelRowMeans}(K)\). Then \(\operatorname{kernelDown}(K)_{ij} = \rho_i\), that is \(\rho \, \mathbf{1}_m^{\top}\). The row means spread across the columns.kernelAcross(k)— Let \(K \in \mathbb{R}^{n \times m}\) and \(\gamma = \operatorname{kernelColMeans}(K)\). Then \(\operatorname{kernelAcross}(K)_{ij} = \gamma_j\), that is \(\mathbf{1}_n \gamma^{\top}\). The column means spread down the rows.
The softmax
Computed naively, exp(1000) overflows to infinity and the softmax becomes NaN. Subtracting a constant
from every score in a row leaves the softmax exactly unchanged, so each row’s maximum is subtracted first;
the log-sum-exp adds it back after the logarithm.
softmaxRows
softmaxRows(z) · z scores n×K → matrix n×K
Let \(Z \in \mathbb{R}^{n \times K}\) (the argument z) and \(m_i = \max_{k} Z_{ik}\). Then
so every row lies in the probability simplex \(\Delta^{K-1}\).
The softmax of each row. Every row is non-negative and sums to one, and scores in the thousands are handled without overflow.
Intuition Turn arbitrary scores into probabilities that keep their order: exponentiate so everything is positive and a gap in score becomes a ratio of odds, then divide by the row total. Only differences between scores matter, which is why subtracting the row maximum first changes nothing but prevents overflow.
Use it when you convert classifier logits into class probabilities, compute attention weights, or define a Boltzmann (softmax) policy over action values — divide the scores by a temperature first to make choices sharper or flatter. Each row is one distribution; for a single score vector pass a 1×K matrix. For the log of the normaliser, use logSumRows rather than ln of a sum.
classificationexplorationnumerical-stability
def softmaxRows(z) ≔ exponentials(z) / lifted(z, einsum("ij->i", exponentials(z)))
logSumRows
logSumRows(z) · z scores n×K → vector n
Let \(Z \in \mathbb{R}^{n \times K}\) (the argument z) and \(m_i = \max_{k} Z_{ik}\). Then
The log-sum-exp of each row, computed stably by taking the row maximum out before the exponentials and adding it back after the logarithm. This is the normaliser a cross-entropy needs.
Intuition A soft maximum: close to the largest score when one dominates, and up to \(\ln K\) above it when all are equal. It is also the log of a sum of probabilities stored as logs, done without leaving log space.
Use it when you compute a softmax cross-entropy (\(\text{logSumRows}(z)_i - z_{i,y_i}\) is the loss for true class \(y_i\)), a log-likelihood that marginalises over hidden states, or a log partition function. Computing ln of a sum of exp overflows at scores near 710; this does not.
numerical-stabilityclassificationprobabilistic-modelling
def logSumRows(z) ≔ ln(einsum("ij->i", exponentials(z))) + rowMost(z)
Building blocks
steadied(z)— Let \(Z \in \mathbb{R}^{n \times K}\). Then \(\operatorname{steadied}(Z)_{ij} = Z_{ij} - \max_{k} Z_{ik}\). Each row with its maximum subtracted.exponentials(z)— Let \(Z \in \mathbb{R}^{n \times K}\). Then \(\operatorname{exponentials}(Z)_{ij} = e^{\operatorname{steadied}(Z)_{ij}} = e^{Z_{ij} - \max_k Z_{ik}}\). The exponential of the shifted scores.
Examples
Matrices are printed row by row.
Distances to centres, and the nearest one
open Prelude.Kernels
fn distances(x, c) ↦ distance(x, c)
fn nearestMask(x, c) ↦ lowest(distance(x, c))
fn nearestDistance(x, c) ↦ markedRow(distance(x, c), lowest(distance(x, c)))
using Examples;
using MLambda.Hilbert.Runtime;
Tensor x = new double[,] { { 0, 0 }, { 3, 4 }, { 1, 1 } };
Tensor c = new double[,] { { 0, 0 }, { 3, 0 } };
var d = Kern.Distances(x, c);
var nearest = Kern.NearestMask(x, c);
var gap = Kern.NearestDistance(x, c);
Output (verified, rounded to 6 significant digits):
d = [[0, 3], [5, 4], [1.41421, 2.23607]]
nearest = [[1, 0], [0, 1], [1, 0]]
gap = [0, 4, 1.41421]
The point (3, 4) is 5 from the origin and 4 from (3, 0), so its mask marks the second centre.
Ties mark every position
open Prelude.Kernels
fn tieMask(d) ↦ lowest(d)
fn farthestMask(d) ↦ highest(d)
fn rowMin(d) ↦ rowLeast(d)
using Examples;
using MLambda.Hilbert.Runtime;
Tensor d = new double[,] { { 2, 1, 1 }, { 0, 5, 3 } };
var low = Kern.TieMask(d);
var high = Kern.FarthestMask(d);
var least = Kern.RowMin(d);
Output (verified):
low = [[0, 1, 1], [1, 0, 0]]
high = [[1, 0, 0], [0, 1, 0]]
least = [1, 0]
Kernel matrices
open Prelude.Kernels
fn gaussian(x, c, g) ↦ rbf(x, c, g)
fn dotKernel(x, c) ↦ linear(x, c)
fn polyKernel(x, c, g, r, d) ↦ polynomial(x, c, g, r, d)
fn laplaceKernel(x, c, g) ↦ laplacian(x, c, g)
using Examples;
using MLambda.Hilbert.Runtime;
Tensor x = new double[,] { { 0, 0 }, { 3, 4 }, { 1, 1 } };
Tensor c = new double[,] { { 0, 0 }, { 3, 0 } };
var k1 = Kern.Gaussian(x, c, (Tensor)0.5);
var k2 = Kern.DotKernel(x, c);
var k3 = Kern.PolyKernel(x, c, (Tensor)1.0, (Tensor)1.0, (Tensor)2.0);
var k4 = Kern.LaplaceKernel(x, c, (Tensor)1.0);
Output (verified, rounded to 6 significant digits):
k1 = [[1, 0.011109], [0.000004, 0.000335], [0.367879, 0.082085]]
k2 = [[0, 0], [0, 9], [0, 3]]
k3 = [[1, 1], [1, 100], [1, 16]]
k4 = [[1, 0.049787], [0.006738, 0.018316], [0.243117, 0.106878]]
For the point (1, 1) against the origin the squared distance is 2, so the RBF value is \(e^{-0.5 \cdot 2} = e^{-1}\); against (3, 0) the inner product is 3, so the polynomial value is \((3 + 1)^2 = 16\).
Softmax on large scores
open Prelude.Kernels
fn probabilities(z) ↦ softmaxRows(z)
fn logNormaliser(z) ↦ logSumRows(z)
using Examples;
using MLambda.Hilbert.Runtime;
Tensor z = new double[,] { { 0, 0, Math.Log(2) }, { 1000, 1000, 1000 } };
var p = Kern.Probabilities(z);
var lse = Kern.LogNormaliser(z);
Output (verified, rounded to 6 significant digits):
p = [[0.25, 0.25, 0.5], [0.333333, 0.333333, 0.333333]]
lse = [1.38629, 1001.10]
The first row’s exponentials are 1, 1 and 2, so the probabilities are quarters and halves and the log-sum-exp is \(\ln 4\). The second row would overflow if computed naively; here it gives \(1000 + \ln 3\).
Centre a kernel matrix
open Prelude.Kernels
fn centred2(k) ↦ centreKernel(k)
using Examples;
using MLambda.Hilbert.Runtime;
Tensor k = new double[,] { { 2, 1 }, { 1, 2 } };
var centred = Kern.Centred2(k);
Output (verified):
[[0.5, -0.5], [-0.5, 0.5]]
Every row mean, column mean and the overall mean are 1.5, so each entry loses 1.5; the rows and columns of the result sum to zero.
A shaped zero
open Prelude.Kernels
fn zerosLike(a) ↦ zeroed(a)
using Examples;
using MLambda.Hilbert.Runtime;
var z = Kern.ZerosLike((Tensor)new double[,] { { 5, 6 } });
Output (verified):
[[0, 0]]