Equilibrium Networks

Imagine that you were a network security researcher working to develop a traffic characterization or anomaly detection algorithm for a packet-level (versus flow-level) security platform. If you were convinced that complex packet interarrival behavior was intrinsic to network traffic, you might not even pursue the use of tools leveraging Poisson behavior. Even armed with the time rescaling theorem, the only reasonable way to estimate packet rates over short time intervals would require using a model like a Poisson cluster process, and you’d still need a good model for the clusters themselves before you felt comfortable. Even though this looks doable, it’d be complicated, and going any further to produce an overall algorithm would just be that much harder. (Although there is a known large-deviation result for Poisson cluster processes, it is considerably more complicated than the result sketched earlier. See Bordenave and Torrisi’s paper “Large deviations of Poisson cluster processes” for more details.) Add to this the fact that most network traffic security platforms don’t normalize their data in a way that is mathematically convenient, and it’s not hard to see why these sorts of techniques have been underutilized.

Other methods besides the one described earlier in this series of posts, including methods based on Markov processes, are also well-suited for providing scalable, long-term solutions to computer network traffic characterization. Equilibrium’s platform provides a solid foundation for incorporating these (and other) techniques. For example, the Girsanov theorem can be applied to produce a simple, scalable likelihood ratio test rooted in martingale theory that can help to determine changes in network traffic patterns with only trivial dependence on any statistical properties of traffic. It can be evaluated on a per-packet basis at high speeds, or periodically at whatever speeds the computational resources available will support.

More complicated large deviation results have also been studied, including one based on the transient fluctuation theorem (TFT) in statistical physics. Unfortunately, the computation of the sample mean entropy production rate that is at the heart of this result is so involved that pursuing a scalable characterization or detection algorithm based on it doesn’t look worthwhile unless the computational requirements can be improved. (See Jiang, Qian, and Qian’s Mathematical Theory of Nonequilibrium Steady States: On the Frontier of Probability and Dynamical Systems for background on this.) Interestingly, however, the TFT turns out to be a corollary of the Girsanov theorem but was discovered decades later. Like the path integral, whose foundations were laid by Wiener in 1923 (when Feynman was five years old), the TFT was effectively rediscovered by physicists long after probability theorists had covered much of the same ground.

The Azuma-Hoeffding inequality provides another powerful tool for analyzing martingales with bounded differences that is similar to large-deviation techniques. Getting good a priori bounds in order to apply this inequality is tricky for martingales derived from Poisson or continuous-time Markov processes, principally because such processes can have a large number of arrivals or transitions in very short time intervals. Even though this “explosion” is very unlikely at any given time, on the relatively rare but inevitable occasions when it happens, it makes a big difference. Though for martingales derived from the Dynkin formula a trick using time-ordered exponentiation can help, there is a high computational cost incurred and it is not clear that the results can be useful without frequent resets anyway. You can cheat and get the bounds after the fact, but you inevitably pay a price for that. (The technical reason for this is that the implied a posteriori measure is not absolutely continuous with respect to the actual measure.) One of our team has actually numerically evaluated the price paid for cheating this way in the case of Poisson processes: the price isn’t always that bad, but sometimes it is, and it’s always there.

Striking the right balance between effectiveness, simplicity, and rigor for building on results like these is something we like to focus on at Equilibrium as part of providing the best possible foundation for a long-term network monitoring solution.

To sketch how the large-deviation approach discussed in the last post in this series can actually work, consider a process with conditional rate function 1 for the first arrivals and thereafter. Take the maximum-likelihood rate estimate over the previous arrivals and see how the resulting putative martingales behave using the following MATLAB code:

function y=blogf(q,N,n,plotflag)

% Piecewise-constant inhomogeneous Poisson
% process at unit rate for N steps, then at
% rate q for N steps. MLE of rate obtained
% at each timestep using a window of n steps,
% then time rescaling theorem applied to
% obtain a putative unit-rate Poisson process.
% Finally the estimated martingale is formed
% and a large-deviation result applied
% to obtain the probability of the observed
% behavior at the beginning and middle
% (i.e., where the jump occurs).

% NB. There is no error checking done (e.g.,
% for n, N) and the code is far from perfect,
% but should be OK for illustrative purposes

% unit rate Poisson process
at1 = cumsum(-log(rand(1,N)));
% rate q  
at2 = cumsum(-log(rand(1,N))/q);    
at = [at1, at1(end)+at2];   % arrival times
iat = [at(1), diff(at)];    % interarrivals

% empirical rate estimate over last n arrivals
qe = zeros(1,2*N);    
for j = 1:n
 qe(j) = j/at(j);
end
for j = (n+1):2*N
 qe(j) = n/(at(j) - at(j-n));
end

% time rescaled process w/ empirical rate
tau = zeros(1,2*N);   
tau(1) = iat(1)*qe(1);
for j = 2:2*N
 tau(j) = tau(j-1) + iat(j)*qe(j);
end

% putative martingale from empirical rescaling
m = (1:2*N) - tau;  

t = 1;  % timescale
% j1 and j2 are respectively the indices
% of the first arrival further than n*t from 0
% and from the Nth arrival
j1 = find(tau > n*t, 1, 'first');
j2 = find(tau > tau(N)+n*t, 1, 'first');
% martingale differences
dm1 = m(j1) - m(1); dm2 = m(j2) - m(N);
% a, x, yy vars from large-deviation result
a1 = dm1/n; a2 = dm2/n;
x1 = a1/t + 1; x2 = a2/t + 1;   
yy = t*n;   % (note: y is output var.)
% log of asymp. prob. of fluctuations
logp1 = yy*(-x1*log(x1) + x1 - 1);   
logp2 = yy*(-x2*log(x2) + x2 - 1);

y = [logp1 logp2];

if plotflag
 figure;plot(tau,m,'k');hold;
 mmmm = [min(m) max(m)];
 plot(tau(N)*ones(1,2),mmmm,'k:');
 axis tight
 xlabel('\tau'),ylabel('N_\tau - \tau')
 s1 = 'actual rate changes from 1 to ';
 s2 = sprintf('%g',q);
 s3 = ' at ';
 s4 = sprintf('%d',N);
 s5 = 'th arrival; MLE rate over ';
 s6 = sprintf('%d',n);
 s7 = ' arrivals';
 title([s1 s2 s3 s4 s5 s6 s7]);
end

Two typical outputs are shown below.

Typical graphical output of blogf for

As above, except .

The following code snippet (which omits the necessary plotting commands) can be used inline for a quick evaluation of the behavior of the overall technique as a function of and

% L = number of runs to average over
N = 10000; L = 40;
q = .65:.05:1.65; n = ceil(logspace(1,3,25));
for j = 1:21
 qj = q(j);
 for k = 1:25
  nk = n(k); temp = zeros(L,2);
  for l = 1:L
   temp(l,:) = blogf(qj,N,nk,0);
  end
  z1(j,k) = mean(temp(:,1));    % initial
  z2(j,k) = mean(temp(:,2));    % middle
 end
end

Using this, it’s easy to see the difference between the log-likelihood of observed fluctuations before and during an abrupt rate change.

Average log-likelihood of the initial steps over 40 observed runs as a function of and . Note that there should be no dependence on in this case, and that in fact no such dependence is evident.

Average log-likelihood of the steps immediately following a rate change over 40 observed runs. Note the tame behavior around .

From the figures, it’s pretty clear that this kind of mathematically oriented technique can work effectively to detect and estimate the likelihood of even relatively small but abrupt rate changes quickly, without relying on heuristics.

But it can do more. Suppose the rate is piecewise linear (for convenience, as a function of arrivals), changing from 1 to linearly over arrivals and then remaining there for another arrivals. The relevant changes to the MATLAB code basically consist of using the lines

qq = linspace(1,q,N);
at2 = cumsum(-log(rand(1,N))./qq);
at3 = cumsum(-log(rand(1,N))/q);
at12end = at1(end)+at2(end);
at = [at1, at1(end)+at2, at12end+at3];

and changing instances of to along with revised plotting commands.

Here the size of the window used makes a big difference:

Typical graphical output of the modified blogf code for .

Typical graphical output of the modified blogf code for .

The average likelihood after the rate begins to change for is approximately for and for . Of course, the rate estimation window needn’t equal the window used for measuring the likelihoods of fluctuations, but this just adds another variable to the analysis.

The overall point is that as arrivals get closer to Poissonian, so that an arrival rate can be well-defined and estimated, the picture for characterization and detection algorithms improves substantially. The next and final post in this series will give a little bit of perspective.