Neural Nets and Anomaly Detection pt 1

More playing with anomaly detection, this time with Keras. All the data is synthetic and can be generated from the code on my GitHub here.

The Data

Our data looks like this:

Typical synthetic data

The top row shows noisy data that have a clear global periodicity. The bottow row shows noisy data with each column having a periodicity independent of the other columns.

Can a Variational Auto Encoder tell the difference?

The Model

The code for our neural net was adapted from this StackOverflow answer here. Very quickly (roughly 10 epochs), the VAE was able to differentiate which group a given square was in with an accuracy of about 90%.

The data projected onto the VAE's bottleneck

Note that I am not going to tune the VAE to get the last drop of accuracy. I'm happy with about 90%.

Note that the trick with VAEs is that the data with which we train it is also the data that we expect. This is simply because we're asking the neural net to reproduce what it's been given. It won't be exact at all since we're deliberately compressing the data through a bottleneck. This is the very definition of an auto encoder. The 'variational' term comes from the fact that the neural net learns the probabilities of the distribution representing the data [Quora].

The Results

The trouble with VAEs is that although you might be able to print a pretty graph of the rendered data, you might need some algorithm to differentiate all the data for you. I'm going to use SciKitLearn to run a KMeans algorithm since I know there are exactly 2 groups (that is, k=2 here).

This is fine in this example where I have two groups of 256 elements. What if we're faced with imbalanced data and are looking for outliers? KMeans does not do so well with that. It doggedly tries to find two clusters of roughly equal size giving an accuracy of about 50% - the monkey score.

Spot the outlier

To this end, we can use the DBScan algorithm, again from SKLearn. This finds my outlier but does take a bit of tuning.

Still playing...

Hypothesis Testing Notes

Recipe

A hypothesis test is only a partial answer to the question "am I being fooled by randomness?". It does not address sampling bias nor measurement errors both of which might be responsible for the apparent effect. [1]

The way to put a hypothesis test together is:

1. Define a test statistic.

2. Define a null hypothesis.

3. Ensure the model can generate synthetic data.

4. Calculate the p-value by generating data (step 3) that fulfills the test statistic.

"All hypothesis tests fit into this framework." [2]

p-value

"Ideally we should compute the probability of seeing the effect (E) under both hypotheses; that is P(E | H0) and P(E | HA).  But formulating HA is not always easy, so in conventional hypothesis testing, we just compute P(E | H0), which is the p-value" [3]

Given our recipe above, this is how we'd calculate the p-value of a (potentially) crooked dice.

1. Run 1000 simulations
2. Find the chi-squared value for each
3. Count the number above a threshold defined by the chi-squared value for the (expected, actual) tuple.
4. Divide by the 1000 trials and you have the p-value

Note that "the p-value depends on the choice of test statistic and the model of the null hypothesis, and sometimes these choices determine whether an effect is statistically significant or not." [4] "If different models yield very different results, that's a useful warning that the results are open to interpretation." [2]

Worked Example

Problem: In ThinkStats2, the assertion is made that first babies arrive earlier than subsequent babies. Given pregnancy duration for 4413 first babies and 4735 subsequent births, there is a 0.078 week different between the two means. Is this significant?

Methodology: we repeatedly sample pairs of births and measure their difference. After N samples, we calculate our mean difference.

With the null hypothesis, we assume there is no effect and sample from all data without distinguishing between the two labels. We do this for N iterations by shuffling the original data and assigning each datum to one of two groups. Then, we take the difference in the means of those two groups. After N samples, we calculate our mean difference.

1 Statistical Inference is only mostly wrong - Downey

2 There is still only one test - Downey

3 There is only one test! - Downey

4 Think Stats2 - Downey

Kleislis

What are they?

Accordingly, the functions that we try to compose are actually A => Monad[B] . A wrapper around them, a category that is naturally associated with a Monad[B], is called a Kleisli category. A => Monad[B] or Kleisli arrows is just a way to compose these sort of functions, nothing more... 

Fun fact: If you compose Kleisli arrows for IO monad, you will get a description of your computer program. Your computer program is essentially one gigantic Kleisli arrow, with some input and output of Unit that acts as a description, and a runtime environment that executes this program works as an interpreter.

[Pavel Zaytsev on Medium.com]

Rob Norris @tpolecat Nov 20 20:39

Kleislis compose, the equivalent bare functions don't. It's just a matter of convenience.
Kleisli[F, A, B] and A => F[B] are isomorphic, pick whichever works better.

John D Cook has a blog post here how the use of Kleislis could have saved the Mars lander from failing (TL;DR; the conversion from kilometers to miles was done twice).

Why are they useful?

They compose whereas monads (usually) don't.

If you know how to map over A[X] and over B[X] you also know how to map over A[B[X]]. Automagically. For free. 

This is untrue for Monad: knowing how to flatMap over A[X] and B[X] doesn’t grant you the power of magically deriving a flatMap for A[B[X]]. 

It turns out this is a well known fact: Monads do not compose, at least not generically.

[Blog of Gabriele Petronella]

From Rob Norris' talk, Functional Programming with Effects

[Rob's talk is here].

Using Haskell notation: "there is a way to compose functions like a -> m b and b -> m c into a -> m c (like the function here, from String to Set[String]). They are called Kleisli functions" [SO]

Ryan Peters [@sloshy on Gitter Jul 14 18:21 2020] describes why he likes Kleisli's for Dependency Injection:

damnit I keep stumbling into Kleisli by accident... 

My architecture as of late has looked something like this:

Config object defining how to read in config arguments/env (using Decline or Ciris)

Dependencies case class that has a constructor which depends on that Config and tries to read it in, converting those values to all of my dependencies.

Anything that needs different dependencies, like concurrent jobs or other processes, is defined as a Kleisli as mentioned.

Main loads a single Resource[IO, Dependencies] and converts all of my Kleisli jobs into Kleisli[IO, Dependencies, ExitCode] and runs it at once.

Kleisli is intimately related with "first do X then do Y" (sequencing) [Gavin Bisesi Gitter May 06 21:02, 2020]. Also, note that ReaderT == Kleisli [Adam Rosien Gitter May 06 19:28 2020]

Usage in HTTP4S

[Gabriel Volpe @gvolpe Jul 06 16:17 2020]

So, HttpRoutes is an alias for 

Kleisli[OptionT[F, *], Request[F], Response[F]]

We say 

def handle(routes: HttpRoutes[F]): HttpRoutes[F] = ...

So that's the type we are expecting to have (we are annotating our function [with] the expected type), meaning we don't need to do the same again within the Kleisli block because the types can be inferred. It would be a different story if we wouldn't have told the compiler what the expected type is.

Then we have the implementation:

Kleisli { req =>

      OptionT {

        A.handleErrorWith(routes.run(req).value)(e =>
          A.map(handler(e))(Option(_)))

      }

    }

Kleisli { ... } is the same as calling Kleisli.apply. So [within] this block, we expect the type to be 

Request[F] => OptionT[F, Response[F]]

Remember that req can be inferred by the compiler, so we know we have req: Request[F]

Then again we call OptionT { ... } which is the same as OptionT.apply. So [within] this block, we expect the type to be 

F[Option[Response]].

Now we get to the final block

A.handleErrorWith(routes.run(req).value)(e => 
   A.map(handler(e))(Option(_))
) 

Let's recap on the types we have

routes:  HttpRoutes[F]

req:     Request[F]

handler: E => F[Response[F]]

A:       ApplicativeError[F, Throwable]

A good exercise would be to split every single part to understand what the types are

routes.run(req) effectively runs the Kleisli by feeding a Request[F], which gives us an OptionT[F, Response[F]]. So we call .value on it to get F[Option[Response[F]]

Now here's the type signature of handleErrorWith

def handleErrorWith[A](fa: F[A])(f: E => F[A]): F[A]

fa is F[Option[Response[F]] for us

so the second part should be 

E => F[Option[Response[F]] 

but notice that handler is defined as 

E => F[Response[F]] 

so we perform that final map(Option(_)) to lift that Response[F] into an Option[Response[F]].

A worked example using HTTP4S

I found this sample on GitHub informative. It presents an HTTP server using Http4s, a database from H2 (from a quick sbt 'show dependencyClasspathFiles' or sbt dependencyTree), Quill and even Flyway to update the database schema. All of this in less than 150 lines!

The whole app is basically this line:

  private val httpApp: Kleisli[IO, Request[IO], Response[IO]] = (videos <+> tags <+> videotags).orNotFound 

where videos, tags and videotags are all HttpRoutes which is just an alias as defined above by Gabriel.

But what is this <+> magic? Well, we're combining Kleisli's by treating their content as semigroups. The results seemed a little counter intuitive to me to begin with. Run this code:

  val add1      = Kleisli[List, Int, Int](x => List(x + 1))

  val divide10  = Kleisli[List, Int, Int](x => List(x / 10))

  val add1Div10 = add1 <+> divide10

  println(add1Div10.run(99))

and you'll see:

List(100, 9)

Straight forward so far. The call to Kleisli.run executed both functions on the input and used an implicit semigroup to glue the answers together within the List.

But HttpRoutes does not deal with Lists, it deals with OptionTs - a monad transformer for Options. How they behave is demonstrated in my GitHub project. Basically, the binary operator defers to OptionTSemigroupK and it looks like this:

def combineK[A](x: OptionT[F, A], y: OptionT[F, A]): OptionT[F, A] = x.orElse(y)

For example, if our F was a List[Int] (that is, we're gluing together OptionT[List, Int]s) then calling <+> on them will result in a Kleisli whose value is a list with just the first Some (or a single None if they're both Nones). 

Note that if our F were IO then the effect associated with any but the first OptionT would not be run (as demonstrated here in my GitHub). In our HTTP app, the flow short circuits on the first combined function that gives a satisfactory answer.

But now we're getting away from Kleisi's and into semigroups. Suffice to say, composing the Kleisli will glue the results together according to the laws of their semigroup.

Bootstrapping

When trying to evaluate which method most increases the accuracy of a neural net, I had previously used a Bayesian method. But, with a technique known as bootstrapping, I could calculate the confidence interval using classical methods.

"The bootstrap is a method for estimating standard errors and computing confidence intervals" [1]. It's very useful when you have small amounts of data. Using it, we resample (with replacement) the original data many, many times until a histogram of our results is more informative.

"It must be noted that increasing the number of resamples, m, will not increase the amount of information in the data. That is, resampling the original set 100,000 times is not more useful than only resampling it 1,000 times. The amount of information within the set is dependent on the sample size, n, which will remain constant throughout each resample. The benefit of more resamples, then, is to derive a better estimate of the sampling distribution." [TowardsDataScience]

Recall that I only had 9 data points for my two methods (L1 and L2 normalisation of the data). But using the bootstrap technique, this lumpy data could be transformed like this for L1 normalization:

Bootstrapped
accuracy using L1 normalization

and this for L2  normalization:

Bootstrapped 
accuracy using L2 normalization

Note how a sparse distribution becomes much more informative.

By excluding the lowest 5% and the highest 5% of results from bootstrapping we can calculate the confidence interval such:

• the 90% confidence interval for L1 normalization was [0.9466, 0.9504].
• the 90% confidence interval for L2 normalization was [0.9412, 0.9456].

[Python code lives in my GitHub]

Note some caveats with confidence intervals:

"... a 95% confidence interval does not indicate that the parameter of interest has a 95% probability of being within the interval.  Ironically, the situation is worse when the sample size is large.  In that case, the CI is usually small, other sources of error dominate, and the CI is less likely to contain the actual value." Statistical Inference is only mostly wrong - Downey

Also note that boostrapping depends on the Central Limit Theorem - the idea that sampling the mean of most distributions eventually forms a Gaussian. This is not universally true. For instance, the Cauchy distribution is of the form (1+x²)^-1 which you might recognise as a standard integral that can be integrated from -∞ to ∞ (giving a function dependent on arctan) so it can be a valid probability distribution. However, if we tried to find the mean by integrating x/(1+x²) (try using integration by substitution and note that it is symmetric around the y-axis) we get an infinity. Hmm.

Indeed, if we change the Python code to generate data from a Cauchy so:

from scipy.stats import cauchy
...
    data = cauchy.rvs(size=n) 

then don't be surprised if you see something like this:

Bootstrapping a Cauchy

Which is not very useful at all.

[1] "All of Statistics" Wasserman.

Notes on Comonads

This excellent post by Eli Jordan can be summarized as comonads offer counit and cojoin, the opposites of on unit and join (or flatten as Cats calls it) monads. 

Why this is useful is that it can model a use case where there is the notion of focus. To demonstrate this, Jordan has created an implementation of Conway's Game of Life (code). Here, the focus is a cell which we determine is either on or off given its neighbours. (Another example of how comonads use focus can be found here where we traverse a graph, node by node).

The case class on which we will call counit and cojoin is an abstraction, Store[S, A], but for the Game of Life example, it is the grid on which the changing cells live. Notably, Store is instantiated with two arguments, 

1. lookup: S => A
2. index: S  

Specifically for our game, these represent:

1. a mapping from a co-ordinate (S) to a Boolean (A) representing whether the corresponding cell is turned on or off.
2. the co-ordinates of the cell that has the focus. That is, S for us is of type (Int, Int).

For each step in our game, we ask our grid (Store) what the neighbourhood is like (lookup) for its focused cell (index). 

We do this by calling its Cats coflatMap function (that we added via a type class) that takes a grid and returns a Boolean. 

This is where it gets interesting. You'll notice that Store has the potential for a partially applied constructor. That is, it looks like:

case class Store[S, A](lookup: S => A)(val index: S) { ...

So, if we instantiate Store with just a S => A, we actually have a S => Store[S, A].

In our Game of Life, this means that for the coflatMap branch of the code, the grid is not actually a mapping to whether the cell is turned on or off (F[Boolean]) but to the grid of the previous iteration on which we will apply Conway's rules to derive Boolean on/off  mappings. That is, it's our cojoin of type of F[F[Boolean]] that needs to be flattened.

Memoization

Recalculating Conway's rules recursively is expensive. Without memoization, you'll see deep stacks. Using jconsole, I saw quite a few:

- eliminated <owner is scalar replaced> (a Store) at Store.counit(Life.scala:69)

which apparently "represents a redundant lock that has been removed in the bytecode." [SO]. This is because the code on each lookup is trying to recalculate the grid from scatch. Memoization prevents this by caching the grid of the last iteration.

What triggers this recalculation is the rendering of the 20x20 grid. This grid obviously leads to 400 evaluations. But the grid (Store) of each iteration has a reference its predecessor's 400 on which it must calculate Conway's rules. Not only do we not want to recalculate that grid from the game's beginning, we don't want to perform a further 3600 (=400 * 8 neighbours) calculation for cells we've already evaluated - hence the caching.

DR in an Azure Cloud Ecosystem

In our topology, we have a Databricks Spark Structured Streaming job reading from an HDInsights Kafka cluster that is VNet injected into our subscription. So, disaster recovery has to take into account two things:

1. the Kafka cluster
2. the blob storage where the data is landed by SSS.

Kafka

This Microsoft document outlines the usual way in securing your data in Kafka (high replication factors for disk writes; high insync replicas for in-memory writes; high acknowledgement factor etc). In addition, an HDInsight cluster are backed by managed disks that provide "three replicas of your data" each witihin their own availability zone that's "equipped with independent power, cooling, and networking" [Azure docs].

So, within a region, things look peachy. Now, how do we get these Kafka messages replicating across region? The HDInisghts documentation suggests using Apache MirrorMaker but note one critical thing it says:

"Mirroring should not be considered as a means to achieve fault-tolerance. The offset to items within a topic are different between the primary and secondary clusters, so clients cannot use the two interchangeably."

This is worrying. Indeed, there is a KIP to make MirrorMaker 2 fix this problem and others like it (like differences in partitions within topics of the same name; messages entering infinite loops etc). Confluent is pushing its Replicator that (it claims) is a more complete solution (there's a Docker trial for it here). And, there is Brooklin, but Azure says this would be self-managed.

Spark and Blobs

At the moment, all the data goes into one region. The business is aware that in the event of, say, a terrorist attack of the data centres, data will be lost. But even the infrastructure guys have the data being replicated from one region to another, note this caveat in the Azure documentation (emphasis mine):

"Because data is written asynchronously from the primary region to the secondary region, there is always a delay before a write to the primary region is copied to the secondary region. If the primary region becomes unavailable, the most recent writes may not yet have been copied to the secondary region."

But let's say that nothing so dramatic happens. Let's assume our data is there once we get our region back on its feet. In the meantime, what has been landed by SSS in the backup region is incompatible with what came before. This is because Spark stores its Kafka offsets in a folder in Gen2. It's just as well Spark is not writing to the directory that the erstwhile live region was using. If we had been writing to a directory that was common to both regions, some finagling would have to be done as we points the Spark job at another directory, effecting the RTO if not RPO.

Aside: A disaster of a different kind

In the old days of on-prem Hadoop clusters, you might see a problem where too many files were created. The consequence would be the Name Node goes down. Sometimes deciphering the Azure documentation is hard but this link  says the "Maximum number of blob containers, blobs, file shares, tables, queues, entities, or messages per storage account" has "No limit". 

Hopefully (caveat: I have not tested this in the wild) this problem has gone away.

Cache vs Delta Cache

Continuing the investigation into a lack of cache coherency in Databricks in my previous post,  I've raised the issue with the people at Azure.

In the document a representative from Microsoft pointed me to, there are two caches at play: Delta Cache and Apache Spark cache. "The Delta cache accelerates data reads by creating copies of remote files in nodes’ local storage" and is enabled by default on certain clusters. In effect, you get it for free when you use Databricks 

The Apache Spark cache is familiar one that is invoked by calling Dataset.cache().

My first complaint is that these two caches behave inconsistently. In the Delta Cache, updates are immediately available. But if we use the Apache Spark cache in conjunction with the Databrick's Delta Lake (not to be confused with the orthogonal Delta Cache), the data is frozen in time.

Now, in the Spark world, Datasets are strictly immutable. But Delta Lake "provides ACID transactions ... on top of your existing data lake and is fully compatible with Apache Spark APIs." Well, once I've invoked .cache(), and update the Delta Lake in the standard way:

    df.write
      .format("delta")
      .mode("overwrite")
      .option("replaceWhere", CONDITION_STRING)
      .save(HDFS_DIRECTORY)

(working example is here) then I cannot see my own updates!

What is, perhaps, even more bizarre can be seen when I try to read the data afresh. A second, totally independent call (but using the same Spark session) to read the data like this:

session.read.format("delta").load(dir)

cannot see the data either! And just to show you that there is nothing up my sleeve, running a completely separate process (code here) to read the data shows that the change has indeed been persisted to HDFS. This appears to be because Spark caches the data plus the Delta Lake meta data too and cares not for Databricks' new semantics.

This brings me to my next complaint - the leaky abstraction. Databricks is trying to retrofit ACID transactions on an API that was built with immutability at its core. The claim that it's "fully compatible with Apache Spark APIs" seems not enitrely true.

I have a meeting scheduled with the MS representative who is currently of the opinion that we should just never call .cache(). On Azure Databrics, this does not seem too bad as the Delta Cache seems pretty fast. It sucks for anybody using Delta Lake on just HDFS.

Summary

If you call .cache() on a Spark Dataset while using Databricks Delta Lake format you will not be able to:

• see updates from any other Spark session
• see your own updates to this data
• see the true data even if you re-read it from the underlying data store

unlesss you unpersist the Dataset.

[Addendum: I spoke to Databricks engineer, Sandeep Katta, on 17 December 2020 and he agreed the documentation is misleading. He says he'll generate a docs PR and submit it to the relevant team to make it clear that one should not use .cache when using Delta Lake]