ML Ops

"ML Ops is a set of practices that aims to deploy and maintain machine learning models in production reliably and efficiently" according to Wikipedia. In particular, my job includes "reproducibility and diagnostics of models and predictions" (see the "Goals" section). Despite repeatability being the bedrock of good science, it's surprising how few data scientists appreciate the importance of a pipeline that can reproduce results on demand.

ML Ops is an especially interesting field because, unlike most areas of software and data engineering, the mapping is often not one-to-one. For instance, in a typical banking app, an input (purchase of a financial instrument) leads to just one output (the client now owns a financial instrument) modulo any cross-cutting concerns such as logging. But in ML models, this is not always the case. Including a socioeconomic feature from the data set may change the output from the model in subtle ways. Or it may not at all. In fact, how can you be sure that you really did include it? Are you sure there wasn't a bug?

There seems to be a dearth of good resources for ML Ops, so here are some observations I've made.

Ontologies

Ontologies seem a really nice idea but I've never seen them work. This is partly due to paradoxes that derive from trying to be all things to all men. True life story: the discharge data for some hospitals was 01/01/1900 for about 100 patients. This lead to cumulative length-of-stays to be negative. "OK," says the ontology team. "It's only 100 patients out of 10 million so let's remove them". But the other 99 columns for these 100 patients were fine. So, during reconciliation with external systems, the QA team had a devil of a job trying to find why their numbers did not add up. They had no knowledge of the missing 100 data points that were absolutely fine other than their discharge dates.

Governance

Weekly town halls that encourage questions and maximise the pairs of eyes looking at the data while the data cleaning takes place. The format of this is very much like Prime Minister's Question Time. In addition, a regularly updated blog and threads that all teams can read aid knowledge transfer. What you really don't want is the data engineers working in isolation and not communicating regularly with downstream clients like the data scientists and QAs.

Test data

Fake data can sometimes be a problem. Say, I create some data but want to bias it. I make all people with a socioeconomic score less than X be in one class rather than the other and run a logistic regression model on it. All goes well. Then, one day, it is decided to have the socioeconomic score one-hot encoded to represent poor and everybody else. The threshold for this partition happens to be X. Then my tests start failing with "PerfectSeparationError: Perfect separation detected, results not available" and I don't know immediately why (this SO reply points out that this causes a model to blow up).

Non-determinism

If there are any functions that have non-deterministic results (for example F.first when you aggregate a groupBy in Spark) then you are immediately at odds with your QA team. For example, a health provider has demographic data including socioeconomic status. But this can change over time. They want to reduce patients events to just patients. This means a groupBy on the patient ID. But which socioeconomic indicator do we use in the aggregation? 

Plenty of Data Sets

When pipelining data, it's better to persist more data sets than fewer. This facilitates QA as otherwise it's hard to understand what functions acted on a row. The disadvantage is that a pipeline may take longer to run since we're constantly persisting files to disk. Still, the greater visibility you get is worth it.

There is no tool as far as I know to see the lineage of functions on a single row.

Repeatability

I recommend having non-production transforms that examine the model's output. For instance, I'm currently working on a work flow that looks for insights in health data. I have a transform that counts the number of statistically significant insights that also have a significant effect (how you calculate these thresholds are up to you). The point being, after ever pull request, the pipeline is run and the number of insights are compared to the previous run. If they change significantly, there's a reasonable chance something nasty was introduced to the code base.

Healthy looking (auto-generated) numbers of insights from my models

Without a decent platform, it's hard to have repeatable builds (see here and here for horror stories). Now, I have my issues with Palantir's Foundry but at least it has repeatability built in. Under the covers, it uses plain, old Git for its SCM capabilities. What's even nicer, models can be built on a 'data' branch that corresponds to the Git branch, so you won't overwrite the data created by another branch when you do a full build. 

Mini Driver

You can easily blow up the driver in Spark by calling collect() on a particularly large dataset. But is there any other way to do so?

This was the problem facing me this week: perfectly sensible looking code that resulted in:

Job aborted due to stage failure: Total size of serialized results of 4607 tasks (1024.1 MiB) is bigger than spark.driver.maxResultSize (1024.0 MiB)

Well, the answer is, yes, innocuous joins can overwhelm the driver if they're broadcast as they are distributed to the executors via the driver.

"This is due to a limitation with Spark’s size estimator. If the estimated size of one of the DataFrames is less than the autoBroadcastJoinThreshold, Spark may use BroadcastHashJoin to perform the join. If the available nodes do not have enough resources to accommodate the broadcast DataFrame, your job fails due to an out of memory error." [Databricks]

The solutions is to call explain() and "review the physical plan. If the broadcast join returns BuildLeft, cache the left side table. If the broadcast join returns BuildRight, cache the right side table."

This can be achieved by calling the Dataset.hint(...) method with the values you can find in JoinStrategyHint.hintAliases.

Generalized Linear Models

GLMs provide a unified framework for modeling data originating from the exponential family of densities which include Gaussian, Binomial, and Poisson, among others. It is defined in [1] as:

"A generalised linear model (or GLM) consists of three components:

1. A random component, specifying the conditional distribution of the response variable Y_i (for the ith of n independently sampled observations).
2. The linear predictor - that is a linear function of regressors:
   η_i = α + β₁X_i1+ β₂X_i2+ ...
3. A smoothe and invertible link function, g, which converts the expectation of the response variable, μ_i ≡ E(Y_i), to the linear predictor" [see 2]. So:
   g(μ_i) = g(E[Y_i])  = η_i

Note that the linear predictor allows for interaction effects, where one variable depends on another and vice versa, or curvilinear effects (ie powers of x terms), etc. Note that the right hand side represents a linear combination of the explanatory variables.

The key to understanding this is by asking: what is an exponential family? The probability mass function would look like this:

f(x|β) = h(x) e^{(T(x)g(β)-A(β))}

where x and β are the data and parameters respectively, and h, g, T, and A are known functions. Now, the thing is that you can shoe-horn a few distributions into this form. Take the binomial distribution you were taught at high school:

f(k, n, p) =  ⁿC_k p^k(1 - p)^n-k

with a bit of basic algebra (try it!), you can get it to look like:

f(x|p) = ⁿC_k e^{(x log[p/(1-p)] – n log(1-p))}

Hey, that's the form of the exponential family! Here, we've set p=k/n and

g = log[p/(1-p)]

Interestingly, point 3 (above) says this function for g equals our linear predictor in 2, above, and you can derive the sigmoid/logit function (try it!). 

Note that this is not why we use the logit function in linear regression. Often, the argument is that it must map the linear predictor that ranges from  ±∞ to [0,1] as we're interested in probabilities. Although logit does this, there are an infinite number of equations that do also. So, here is a link that details why only logit can be the only suitable function.

Now, remember that a Bernoulli distribution is a just a binomial when n=1, so this completely describes the case for logistic regression. This is why when we're looking at a model with a binary response, we tell the GLM to use the Binomal family (see, for instace, here in the Spark docs). What we do for other use cases I shall deal with in another post.

[1] Applied Regression Analysis & Generalized Linear Models, Fox, 3rd edition.

Spark Log Likelihood

Log-likelihood

I wanted to see how the number of iterations in a Generealized Linear Model effected the performance of my mode. One way to do this is with log-likelihood. "Probability is used before data are available to describe plausibility of a future outcome, given a value for the parameter. Likelihood is used after data are available to describe plausibility of a parameter value." [Wikipedia]. In Bayesian inference, if P(C|D) is the probability, P(D|C) is the likelihood. We take the logs to make the maths easier.

Spark

My model was built using Spark's ML library. With it, I was trying to model the number of hours spent waiting in A&E. If the patient waits more than 4 hours, they have breached our target. The code looks like this:

from pyspark.ml.regression import GeneralizedLinearRegression
glr = GeneralizedLinearRegression(family="binomial", link="logit", featuresCol=OUTPUT_COL, labelCol=self.dependent_column, predictionCol=PREDICTION_COLUMN)

Since I am modeling the number of hours waiting, a binomial is an appropriate distribution.

Unfortunately, Spark does not appear to have a log-likelihood function out-of-the-box for GLMs. So, I tried to write my own.

In the Python library, StatsModels, the (abbreviated) code for Logit log likelihood is:

def loglike(self, params):
    q = 2*self.endog - 1
    X = self.exog
    return np.sum(np.log(self.cdf(q*np.dot(X,params))))

def cdf(self, X):
    X = np.asarray(X)
    return 1/(1+np.exp(-X))

from discrete_model.py.

This I turned into PySpark code:    

from scipy import special

    def log_likelihood(self) -> float:
        '''
        PySpark version of the code in StatsModel's discreet_model.Probit.loglikelihood function.
        :return: the log likelihood of the data on which this object was fitted.
        '''

        log_likelihood = "log_likelihood"
        inverse = "inverse"
        lin_pred = "lin_pred"
        df = self.df
        params = self.model.coefficients
        intercept = float(self.model.intercept)

        @F.udf
        def dot_params(exog) -> float:
            vecs = exog.toArray()
            assert len(vecs) == len(params)
            x_dot_params = sum(map(lambda x: x[0] * x[1], zip(vecs, params)))
            return float(x_dot_params) + intercept

        @F.udf
        def loglike_obs(endog, mu) -> float:
            n = 1.
            y = endog * n
            x = (special.gammaln(n + 1) - special.gammaln(y + 1) -
                    special.gammaln(n - y + 1) +
                    y * math.log(mu / (1 - mu + 1e-20)) +
                    n * math.log(1 - mu + 1e-20))
            return float(x)

        df = df.withColumn(lin_pred, dot_params(F.col(OUTPUT_COL)))  # + offset_exposure = 0. by default
        df = df.withColumn(inverse, 1. / (1 + F.exp(-1 * F.col(lin_pred))))
        df = df.withColumn(log_likelihood, loglike_obs(F.col(self.dependent_column), F.col(inverse)))
        df = df.agg(F.sum(log_likelihood).alias(log_likelihood))
        rows = df.select(log_likelihood).collect()
        return rows[0][0]

Messy but it works.

AIC

As it happens, I didn't need to bother. Log-likelihood can be derived from a method that is already in the Spark API.

Once you call GeneralizedLinearRegression.fit, you'll have a GeneralizedLinearRegressionModel. An attribute of this is the summary (of type GeneralizedLinearRegressionTrainingSummary). From this you can get the AIC:

>>> model.summary.aic
116384.05296519656

The AIC is the Akaike Information Criteria, a metric that uses Kullback-Leibler divergence. The equation is a function of the log-likelihood and the number of parameters in the model (including the intercept).

Note that the AIC is calculated based on the choice of family when we create our GeneralizedLinearRegression object.

To get the log-likelihood, some code like this:

        aic = model.summary.aic
        k = (data_frame.count() - model.summary.degreesOfFreedom)
        log_likelihood = ((2 * k) - aic) / 2

should suffice.

I'd previously used log-likelihood as a metric with which I could compare the efficacy of "number of iterations" values. Here, the data was the same and so was the model. AIC is useful if we're comparing two (or more) different models.

Impure Evil

"Debugging is like being the detective in a crime movie where you are also the murderer."

Here's a case study of the dangers of using functions that are not pure in a lazily evaluation environment. While it's obvious that the rand function is not pure, it's less obvious with others.

The requirements

We have a table of national accident & emergency events that also contains patient data. We want just the patients from this data. Some patients appear more than once. The question is: how do we reduce multiple rows for the same patient into a single row?

This is not as straightforward as it sounds. The same patient may appear in the events that are years apart. In that time, they may have live in a different area code. So, one solution was to just choose any as it was not overly important. Therefore, we used the first function. After all, why not? We just want an arbitrary value in the event we have more than only one choice.

Once we have this data, we want to express the ratio of patients with property Z to those without segmented on hospital attended. It is not relevent what Z actually is but we will have to make two calculations: one for the set of those with Z and one for all patients irrespective of Z.

The problem

When we ran this in Spark and saw some values where the numerator was greater than the denominator. 

The diagnosis

The problem arises when we use impure functions, especially with lazy evaluation.

1. The first function will return a non-deterministic value for patients with multiple attendances. There is no guarantee which it will choose.
2. Spark is lazy by default.
3. So, upon calculating the denominators, first is invoked and it gives X for a given patient. Then, when calculating the numerators, first is invoked again (because of #2) but this time returns Y (because of #1) for the same patient. The result is that sometimes more people will fit the criteria for the numerator than the denominator for any given slice-and-dice of the data.

Note that this would also be true of most database implementations. For instance, SQL Server's FIRST_VALUE is also non-deterministic.

Solution

Calling cache() will sometimes help but is not robust. If, for instance, an executor crashes, the data is regenerated by another instance and it may not necessarily be the same as what was generated before.

A better solution is to adopt a deterministic algorithm when choosing from a collection of values. We used max but it could be anything else. 

However, note that there is no guarantee that there is a consistency between values. For instance, we have a patient's area code expressed as a code recognised nationally and also as a code recognised by the health providers. There should be a one-to-one mapping between them. But given more than one of each, applying a lexical max to both may not produce a valid pair.

Big Data or Pokeman?

It's hard to keep track of developments in big data. So much so that there is a quiz to see if you can differentiate big data tools from Pokemon characters. It's surprisingly difficult.

So, to make it look like an old fool like me can avoid being dissed by the cool kids, this is a cheat sheet of what's currently awesome:

Trino
Once known as PrestoSQL, Trino is a SQL engine that can sit on top of heterogenous data sources. PrestoSQL is offered by AWS under the name "Athena". The incubating Apache Kyuubi appears to be similar but tailored for Spark.

Amundsen
Amundsen is a "data discovery and metadata engine". Apache Atlas is an Hadoop based metadata and governance application written in Java.

Apache Griffin
Griffin is a JVM-based tool for checking for data quality. TensorFlow Data Validation and Great Expectations for Python.

Apache Arrow
Arrow is a language-agnostic columnar processing framework. You might need it if you want to use User Defined Aggregate Functions in PySpark [StackOverflow]. It's written in a number of languages, predominantly C++ and Java and can help leverage GPUs.

Azkaban
Azkaban is a Hadoop workflow management tool from LinkedIn. It's open source and written in Java.

Koalas
Koalas brings the Pandas API to PySpark. It depends on Arrow to do this, apparently.

DBT
DBT is an open source Python project that does the T in ELT. Transforms are in templated SQL, apparently.

Apache Pinot
Pinot is a distributed, columnar data store written in Java that ingests batched and streaming data. It's a little like Apache Druid. "The only sustainable difference between Druid and Pinot is that Pinot depends on Helix framework and going to continue to depend on ZooKeeper, while Druid could move away from the dependency on ZooKeeper. On the other hand, Druid installations are going to continue to depend on the presence of some SQL database." [Medium]

Debezium
From the docs: "Debezium is an open source [Java] project that provides a low latency data streaming platform for change data capture (CDC). You setup and configure Debezium to monitor your databases, and then your applications consume events for each row-level change made to the database."

Samza

Apache Samza allows real-time analysis of streams. It comes from the same people who gave us Kafka (LinkedIn) and is written in Java and some Scala.

Pulsar

Apache Pulsar is like Apache Kafka but adds messaging on top of its streaming offering. It's mostly written in Java.

AirByte

This appears to be a half Java, half Python collection of connectors to "sync data from applications, APIs & databases to warehouses, lakes & other destinations."

Hudi

A Java/Scala Apache project that  stands for "Hadoop Upserts Deletes and Incrementals" which pretty much describes it. It can integrate with Spark.

DataHub

This is an open source "data catalog built to enable end-to-end data discovery, data observability, and data governance" [homepage] that came out of LinkedIn. It is written in Python and Java.

The adoption within the industry for big data tools can be found in this report [LinkedIn]

Andreesen-Horowitz make their predictions here.

More Logistic Regression

Log-Likelihood
How good is our model? I'm told the best way to think about the likelihood in a logistic regression is as the probability of having the data given the parameters. Compare this to what logistic regression is actually doing: giving the most probable parameters given the data.

First, some terminology:

Exo and Endo
"If an independent variable is correlated with the error term, we can use the independent variable to predict the error term, which violates the notion that the error term represents unpredictable random error... This assumption is referred to as exogeneity. Conversely, when this type of correlation exists, which violates the assumption, there is endogeneity." - Regression Analysis, Jim Frost.

"An endogenous variable is a variable in a statistical model that's changed or determined by its relationship with other variables within the model. In other words, an endogenous variable is synonymous with a dependent variable, meaning it correlates with other factors within the system being studied. Therefore, its values may be determined by other variables. Endogenous variables are the opposite of exogenous variables, which are independent variables or outside forces. Exogenous variables can have an impact on endogenous factors, however." [Investopedia]

Spark
Unfortunately, Spark does not seem to have a calculation for log-likelihood out of the box. This forced me to code my own. I looked at the code in the Python library, StatsModels, and converted it to PySpark.

Comparing the output from Spark was very data dependent. I guess this is inevitable since Spark uses IRLS [Wikipedia] as the solver and StatsModels was using l-BFGS. I was forced to use l-BFGS in StatsModels as I was otherwise "singular matrix" errors [SO].

But then, a passing data scientist helpfully pointed out that I was looking at the wrong StatsModel class (Logit in discreet_model.py not GLM in generalized_linear_model.py). Taking the GLM code, I got within 6% of StatsModels when both run on the same data. Could I do better? Well, first I rolled back the use of l-BFGS to the default solver (which is IRLS for both libraries) and removed a regularization parameter in the Spark code that had crept in from a copy-and-paste - oops. Now, the difference between the two was an impressive 4.5075472233299e-06. Banzai!

We may choose to have no regularization if we're looking for inferential statistics (which I am). This might lead to overfitting but that's OK. "Overfitting is predominantly an issue when building predictive models in which the goal is application to data not used to build the model itself... As far as other applications that are not predictive, overfitting is more secondary" [SO]

The Data
The size of data makes a difference. The more data, the lower the likelihood it is explained by the parameters. This is simply because the state space is bigger and that means there are more potentially wrong states.

Also, adding a bias of a single feature can change the log-likelihood significantly. Making a previously unimportant feature biased to a certain outcome 2/3 of the time reduced the LL in my data by an order of magnitude. Not surprisingly, if the data is constrained to a manifold, then it's more likely your model will find it.