SPINE — User Guide
Scale Partitioning by Incremental Nested Entropy.
SPINE takes a nonnegative measure over a collection of objects and cuts it into scale strata. It finds where the collection separates into distinct scales, and how strongly, without being told how many groups to look for.
You supply two things: a measure, one nonnegative number per object, and an orientation, saying which end of that measure is the atypical one. Those are the scientific input. Everything after them is computed.
Run this guide section by section. Sections 1 to 5 are the basics. Sections 6 onward are worked scenarios, each on a different kind of data, and they can be read in any order.
Author Abd AlRahman AlMomani Assistant Professor of Data Science and Mathematics, ERAU, AZ almomaniar@gmail.com
1. Before you start
Nothing to install. MATLAB's current folder needs to be the MATLAB folder of the SPINE package, the one holding spine.m and this guide. Everything else follows from there.
if isempty(which('spine'))
error(['Change MATLAB''s current folder to the folder containing ' ...
'spine.m, then run this guide again.']);
end
2. Your first call
Two hundred objects drawn from a standard normal cloud, plus twelve sitting a long way out. Nothing tells SPINE that there are two groups.
rng(11);
X = [randn(200,2); randn(12,2) + 9];
[scores, Result] = spine(X, 'dist2mean', "upper");
fprintf('strata found %d\n', Result.nStrata);
fprintf('cut at measure %.4g\n', Result.boundaryValues);
fprintf('objects flagged %d of %d\n', sum(scores > 0), numel(scores));
3. Reading the result
Scores run from 0 for the typical stratum to 1 for the most atypical, one per object, in your object order. Result carries the rest, also in your order.
disp(Result)
% boundaryValues is the result you can carry elsewhere. It is a value in your
% own measure's units, so the same cut can be applied to new data or quoted in
% a paper. Info.boundaries, by contrast, are ranks into an internally sorted
% copy and mean nothing outside this call.
spineSummary says the same partition in words. It recomputes nothing, so it also works on a Result loaded back from a file long after the run.
spineSummary(Result)
4. Seeing it
Two figures, each with its own job.
spineDataPlot draws your data in the space it lives in, coloured by score. Objects are drawn in increasing order of score, so the atypical ones are never buried under the bulk, and the colour bar is ticked at the strata with their populations, so it doubles as the legend.
spineDataPlot(Result, 'Title', 'a bulk and a far group');
spinePlot shows where the cuts came from: the sorted measure with the cut marked, and the incremental nested entropy profile behind it.
spinePlot(Result, 'Title', 'a bulk and a far group');
5. The measure is the modelling step
This is the part worth slowing down for. The same data under two measures can give opposite answers, and both are correct — they answer different questions.
Three groups: a bulk, a halo displaced a little and spread widely, and an island displaced a long way and packed tight.
rng(4242);
bulk = randn(200,2);
halo = 5.0*[cos(2*pi*rand(30,1)), sin(2*pi*rand(30,1))] + 2.0*randn(30,2);
island = [12 12] + 0.25*randn(20,2);
D = [bulk; halo; island];
group = [ones(200,1); 2*ones(30,1); 3*ones(20,1)];
% Run the same engine under four measures and compare what each one calls
% atypical. Nothing changes except the modelling step.
measures = ["dist2mean", "robustrownorm", "knn", "dist2median"];
fprintf('%-16s %7s %8s %8s %8s\n', 'measure', 'strata', 'bulk', 'halo', 'island');
Results = cell(numel(measures), 1);
for j = 1:numel(measures)
[s, R] = spine(D, measures(j), "upper", 'K', 5);
Results{j} = R;
fprintf('%-16s %7d %8.3f %8.3f %8.3f\n', measures(j), R.nStrata, ...
mean(s(group==1)), mean(s(group==2)), mean(s(group==3)));
end
% Read the island column. dist2mean calls it maximally atypical, 1.000. knn
% calls it perfectly typical, 0.000, and picks out the halo instead: the island
% is far away but densely packed, so its nearest-neighbour distances are small.
%
% Distance from a centre and local density are different questions, and the
% measure is where you say which one you meant. Neither answer is wrong.
%
% docs/measure-guide.html is the fuller treatment, and it is the thing to read
% before applying SPINE to your own problem.
Side by side, the difference is visible rather than tabulated.
spineDataPlot(Results{1}, 'Title', 'dist2mean: far from the centre');
spineDataPlot(Results{3}, 'Title', 'knn: in a sparse neighbourhood');
spineCompare packages this whole comparison: one call runs every measure, reports what each found, and prints the pairwise agreement on which objects were flagged, not merely how many.
spineCompare(D, measures, "upper", 'K', 5)
6. Scenario: rare categories
Data structure: a plain vector of counts.
Forty tokens seen once or twice, six seen dozens of times. Here the RARE ones are the interesting ones, so small is atypical and the orientation is "lower".
rng(20);
counts = [randi([1 3], 40, 1); randi([55 70], 6, 1)];
[sFreq, ResFreq] = spine(counts, 'frequency', "lower");
fprintf('strata %d, flagged %d of %d as rare\n', ...
ResFreq.nStrata, sum(sFreq > 0), numel(sFreq));
fprintf('cut at frequency %.4g\n', ResFreq.boundaryValues);
% Getting the orientation backwards does not produce an error. It produces a
% confident answer to a question you did not mean to ask.
6b. Scenario: rare states in symbolic data
Data structure: text. A string array, a categorical array, a cell array of character vectors, or a table mixing text and numbers.
Use 'frequency' when you already hold one count per category. Use 'statefrequency' when you hold raw observations that still need tallying: each ROW is one observed state, rows are matched exactly, and the objects become the distinct states.
Four thousand draws from a sixty-word vocabulary with a heavy head.
rng(9);
nTok = 60;
weights = 1 ./ ((1:nTok).^1.1);
draws = randsample(nTok, 4000, true, weights/sum(weights));
words = compose('w%02d', draws);
[sWord, ResWord] = spine(words, 'statefrequency', "lower");
fprintf('%d observations over %d distinct states\n', ...
numel(draws), numel(sWord));
fprintf('strata %d, states flagged as rare %d\n', ResWord.nStrata, sum(sWord > 0));
Text and numbers can be mixed, and the state is the whole row.
Log = table(["ok";"ok";"warn";"ok";"fault";"warn";"ok";"fault"], ...
[ 0 ; 0 ; 1 ; 0 ; 2 ; 1 ; 0 ; 3 ], ...
'VariableNames', {'status','code'});
[~, ResLog] = spine(Log, 'statefrequency', "lower");
disp([ResLog.objectIds, string(ResLog.measure.counts)]);
% Two cautions. Numeric columns are matched at full precision, so round or bin
% first if near-misses should count as the same state. And the object count is
% the number of DISTINCT states, not of observations: forty thousand rows over
% three states give SPINE three objects, which is not enough to partition.
%
% To score the original rows rather than the states:
%
% rowScores = sWord(ResWord.measure.stateOfRow);
7. Scenario: hubs in a network
Data structure: a square adjacency matrix.
The objects are nodes, the measure is the degree, and the orientation is "upper" because hubs are the atypical ones. Three networks, each asking a different question of the same measure. The generators are local functions at the end of this guide.
A. A small-world backbone with five designated hubs. Tight, regular degrees plus a handful of deliberate hubs, so the ground truth is unambiguous.
W = smallWorld(200, 6, 0.10, 7);
[A1, hubIds1] = addHubs(W, 5, 0.18, 7);
truth1 = false(size(A1,1),1); truth1(hubIds1) = true;
[s1, R1] = spine(A1, 'degree', "upper");
deg1 = sum(A1 ~= 0, 2);
fprintf('A. %d nodes, degrees %d to %d, cut at %g\n', ...
size(A1,1), min(deg1), max(deg1), R1.boundaryValues);
fprintf(' designated hubs found %d of %d, other nodes flagged %d\n', ...
sum(s1 > 0 & truth1), sum(truth1), sum(s1 > 0 & ~truth1));
B. A scale-free network with six designated hubs. Preferential attachment makes hubs of its own accord, so the designated ones are not the only real hubs present.
S = scaleFree(180, 2, 20260810);
[A2, hubIds2] = addHubs(S, 6, 0.22, 20260810);
truth2 = false(size(A2,1),1); truth2(hubIds2) = true;
[s2, R2] = spine(A2, 'degree', "upper");
deg2 = sum(A2 ~= 0, 2);
extra = find(s2 > 0 & ~truth2);
fprintf('B. designated hubs found %d of %d, other nodes flagged %d\n', ...
sum(s2 > 0 & truth2), sum(truth2), numel(extra));
fprintf(' those extra nodes have degrees %s, against a bulk maximum of %d\n', ...
mat2str(sort(deg2(extra)).'), max(deg2(s2 == 0)));
% They are not false positives. Preferential attachment produced genuine hubs,
% and the degree measure does not care which were designed and which grew.
C. The same backbone with no hubs at all. One stratum is the answer: these degrees sit at a single scale. That is a property of the network, not of the cut.
A3 = smallWorld(200, 6, 0.10, 7);
[s3, R3] = spine(A3, 'degree', "upper");
deg3 = sum(A3 ~= 0, 2);
fprintf('C. degrees %d to %d, strata %d, flagged %d\n', ...
min(deg3), max(deg3), R3.nStrata, sum(s3 > 0));
spineNetworkPlot draws the graph itself: a force layout with nodes coloured by score and sized by degree, so hubs read as hubs.
spineNetworkPlot(R1, A1, 'Title', 'small-world with designated hubs', 'Seed', 3);
spineNetworkPlot(R2, A2, 'Title', 'scale-free', 'Seed', 3);
8. Scenario: how many components does this matrix have?
Data structure: a rectangular data matrix, treated spectrally.
Three strong directions buried in noise. The measure is the singular values themselves, so the OBJECTS are the singular values, not the rows.
rng(77);
U = orth(randn(60,60)); Vt = orth(randn(25,25));
sv = [40 26 15, 0.6*ones(1,22)];
M = U(:,1:25) * diag(sv) * Vt.';
[sSv, ResSv] = spine(M, 'singular', "upper");
fprintf('objects (singular values) %d\n', numel(sSv));
fprintf('strata %d, components above the noise floor: %d\n', ...
ResSv.nStrata, sum(sSv > 0));
% Because the objects are no longer the rows, there is no per-object data to
% draw, and spineDataPlot would tell you so. spinePlot still works, and the
% measure panel is exactly the spectrum with the cut on it.
spinePlot(ResSv, 'Title', 'singular values, with the cut');
9. Scenario: a table with identifiers
Data structure: a table with row names.
Fifty-five ordinary sites and five anomalous ones. The identifiers survive the round trip, so you get back the names of what was flagged, not row numbers.
rng(8);
sites = [randn(55,3); randn(5,3) + 7];
T = array2table(sites, 'VariableNames', {'x','y','z'}, ...
'RowNames', compose("site%02d", (1:60).'));
[sT, ResT] = spine(T, 'robustrownorm', "upper");
fprintf('table: %d objects, %d strata\n', numel(ResT.mu), ResT.nStrata);
disp(ResT.objectIds(sT > 0).');
10. Scenario: a struct array of records
Data structure: a struct array, one element per object.
Numeric scalar fields become the features. A text field becomes the identifiers, automatically.
rng(15);
S = struct('label', {}, 'width', {}, 'height', {}, 'mass', {});
for i = 1:45
far = i > 41;
S(i).label = sprintf('part-%03d', i);
S(i).width = 10 + randn + 14*far;
S(i).height = 5 + randn + 9*far;
S(i).mass = 2 + randn + 6*far;
end
[sS, ResS] = spine(S, 'robustrownorm', "upper");
fprintf('struct array: %d objects, %d strata\n', numel(ResS.mu), ResS.nStrata);
disp(ResS.objectIds(sS > 0).');
11. Scenario: three dimensions
Data structure: an n-by-3 matrix.
A dense core, a thin shell around it, and far outliers. Drag the figure to rotate it.
rng(90210);
dirs = @(m) randn(m,3) ./ sqrt(sum(randn(m,3).^2, 2));
V1 = randn(90,3); U1 = V1 ./ sqrt(sum(V1.^2,2));
V2 = randn(18,3); U2 = V2 ./ sqrt(sum(V2.^2,2));
X3 = [0.9*randn(260,3); U1.*(7 + 0.25*randn(90,1)); U2.*(34 + 3*randn(18,1))];
[~, Res3] = spine(X3, 'dist2median', "upper");
fprintf('3-D: %d strata, cut radii %s\n', Res3.nStrata, ...
join(compose("%.3g", Res3.boundaryValues(:).'), ", "));
% The three strata score 0, 0.18 and 1, so on a score colour scale the core and
% the shell come out nearly the same shade. Colouring by stratum spaces them
% evenly instead, and the bar is labelled "stratum" so nothing is misread.
spineDataPlot(Res3, 'Title', 'core, shell, outliers', 'ColourBy', "stratum");
12. Scenario: many features at once
Data structure: a wide matrix, more features than can be drawn.
Twelve features. The measure reduces them to one number as usual; only the drawing is affected, and it shows the first two principal components with the axes labelled to say so.
rng(64);
X12 = [randn(240,12); randn(18,12) + 5];
[~, Res12] = spine(X12, 'robustrownorm', "upper");
fprintf('12 features: %d strata, %d flagged\n', ...
Res12.nStrata, sum(Res12.scores > 0));
spineDataPlot(Res12, 'Title', '12 features, projected');
13. Scenario: you already have the measure
Data structure: a vector that is already the quantity you care about.
If your domain hands you the number — an energy, a claim size, a residual — use it directly with 'asis'. A domain quantity almost always beats a constructed one, because it carries meaning your reader already understands.
rng(5);
energy = [gamrnd(2, 1.0, 300, 1); gamrnd(2, 25, 9, 1)];
[sE, ResE] = spine(energy, 'asis', "upper");
fprintf('asis: %d strata, cut at %.4g, %d flagged\n', ...
ResE.nStrata, ResE.boundaryValues, sum(sE > 0));
14. Scenario: atypical at both ends
Sometimes both extremes matter and the typical objects are in the middle. Bilateral orientation needs an anchor: an object you consider typical.
Here the measure is distance from the median, and the anchor is the object closest to it.
rng(42);
v = [randn(200,1)*1 + 20; randn(8,1)*0.5 + 3; randn(8,1)*0.5 + 46];
dev = abs(v - median(v));
[~, anchor] = min(dev);
[sB, ResB] = spine(dev, 'asis', "bilateral", 'BilateralAnchor', anchor);
fprintf('bilateral: %d strata, %d flagged at the two ends\n', ...
ResB.nStrata, sum(sB > 0));
fprintf('anchor object %d has score %.3f\n', anchor, sB(anchor));
15. Scenario: robust or not
When you are hunting for objects at a different scale, do not standardize by statistics those same objects distort. The standard deviation is inflated by the very group you are looking for; the MAD is not.
rng(20260810);
Xr = [randn(240,8); randn(45,8) + 6; randn(15,8) + 30];
[~, Rz] = spine(Xr, 'zrownorm', "upper");
[~, Rmad] = spine(Xr, 'robustrownorm', "upper");
fprintf('zrownorm %d strata\n', Rz.nStrata);
fprintf('robustrownorm %d strata (the planted structure)\n', Rmad.nStrata);
% Prefer robustrownorm over zrownorm, and dist2median over dist2mean, unless
% you have a reason not to.
16. Scenario: nothing to find
Structureless data gives one stratum, and that is an answer rather than a failure: this measure does not separate these objects into scales. It is a property of the measure you chose, not of the cut.
rng(3);
[~, ResFlat] = spine(randn(300,4), 'dist2median', "upper");
fprintf('unstructured data: %d stratum\n', ResFlat.nStrata);
% The right response is to ask whether a different measure would express the
% structure better, not to loosen anything — there is nothing to loosen.
17. Options you may want
Every option not listed below is passed straight through to spineOptions.
'Objects', "columns" your matrix has objects in columns, not rows
'K', 8 neighbour count for the knn measure
'ObjectIds', ids supply identifiers yourself
'BoundaryMode', "selected" keep a selected subset of the candidate cuts
'Missing', "omit" drop objects with missing values; their scores return NaN
'Thin', true very large collections: locate the cut on a thinned copy
'BilateralAnchor', id the typical object, for bilateral orientation
'ProfileWeights', w weight the profile ranks, selected mode only
An example of BoundaryMode: candidate mode keeps every candidate boundary, selected mode chooses among them.
rng(99);
Xs = [randn(150,2); randn(30,2)+4; randn(10,2)+14];
[~, Rc] = spine(Xs, 'dist2median', "upper");
[~, Rs] = spine(Xs, 'dist2median', "upper", 'BoundaryMode', "selected");
fprintf('candidate mode %d strata, selected mode %d strata\n', ...
Rc.nStrata, Rs.nStrata);
18. Where to go next
../docs/measure-guide.html how to choose a measure, in full
../README.md the reference for everything Result gives you
../paper/figures/ one script per manuscript figure
help spineMeasure lists the ten measures. help spine documents every option. help spineDataPlot, help spineNetworkPlot and help spinePlot cover the figures.
close all;
Local functions
The network generators used in section 7. Both are plain MATLAB and need no toolboxes.
function A = scaleFree(n, m, seed)
% Barabasi-Albert preferential attachment: each new node attaches to m existing
% nodes chosen with probability proportional to their current degree.
rng(seed);
m0 = m + 1;
A = zeros(n);
A(1:m0, 1:m0) = 1 - eye(m0);
pool = repelem((1:m0).', m0-1);
for v = m0+1:n
chosen = zeros(1,0);
while numel(chosen) < m
t = pool(randi(numel(pool)));
if ~any(chosen == t)
chosen(end+1) = t; %#ok<AGROW>
end
end
A(v, chosen) = 1;
A(chosen, v) = 1;
pool = [pool; chosen(:); repmat(v, m, 1)]; %#ok<AGROW>
end
end
function A = smallWorld(n, k, beta, seed)
% Watts-Strogatz: a ring lattice in which each edge is rewired at random with
% probability beta. Degrees stay tightly concentrated around k.
rng(seed);
A = zeros(n);
for i = 1:n
for j = 1:k/2
t = mod(i+j-1, n) + 1;
A(i,t) = 1;
A(t,i) = 1;
end
end
for i = 1:n
for j = 1:k/2
t = mod(i+j-1, n) + 1;
if A(i,t) == 1 && rand < beta
cand = find(A(i,:) == 0);
cand(cand == i) = [];
if ~isempty(cand)
A(i,t) = 0;
A(t,i) = 0;
nt = cand(randi(numel(cand)));
A(i,nt) = 1;
A(nt,i) = 1;
end
end
end
end
end
function [A, hubIds] = addHubs(A, nHub, frac, seed)
% Attach nHub new nodes, each wired to a fraction frac of the existing ones.
rng(seed);
nBase = size(A,1);
A(nBase+nHub, nBase+nHub) = 0;
hubIds = (nBase+1 : nBase+nHub).';
for h = 1:nHub
v = hubIds(h);
targets = randperm(nBase, round(frac*nBase));
A(v, targets) = 1;
A(targets, v) = 1;
end
end