The decomposition of yield curves using principle component analysis (PCA) is a common procedure, allowing for yield curve dynamics to be explained in a small number of underlying drivers. The first three principle components are commonly used and are given economic interpretations as in Litterman and Scheinkman, 1991 and Kavir Patel, 2018.

PC1 - Level: The first principle component tends to have equal values across all maturities, so controls the overall level of the curve.

PC2 - Slope: The second component usually has a positive relationship with maturity, which allows this component to determine the slope of the curve and the bulk of the differences in yield between short and long maturities.

PC3 - Curvature: The third and usually last component tends to determine how bowed the curve is, leading to it being referred to as curvature.

These three components can be visualised by plotting the loadings on each principle component for each maturity, the eigenvectors of the covariance matrix, against maturity, as in Figure 2.

The use of PCA to decompose the yield curve into factors with distinct economic interpretations raises some questions on how the reliability of the decomposition, as the decomposition is com- pletely informed by historical data.

This question is explored by Rodger Lord, 2007, where the authors investigate the circumstances under which the economic interpretations are consistent with their Level, Slope and Curvature des- ignations, concluding that there are some years in which there are no economic interpretations for the components. While Rodger Lord, 2007 use the shape and sign changes in the eigenvectors to determine if the economic interpretation is valid, Hou, 2022 uses a comparison against the model developed by Nelson and Siegel, 1987, which explicitly models the components, to test if the PCA based components are poorly estimated, also concluding that there are time periods where the interpretations of the components are inconsistent with their usual interpretations.

This investigation aims to determine impact of the choice of covariance estimator on the estimation of the principle components of bond yields and whether specific estimators provide more consistency on historic data.

Approach

Two approaches are explored to quantify how well estimated the principle components are:

Sign change of eigenvectors: This borrows the method and sign change rules used by Rodger Lord, 2007, which checks if ”first three factors have respectively zero, one and two sign changes”, in which case the level, slope and curvature are assumed to be estimated in accor- dance with their economic interpretations. Applying these rules to the results in Figure 2 we can see that the decomposition is valid

Visual Inspection: Inspection of graphs of eigenvectors against maturity to try find any insights that may have been missed by either of the other methods.

These approaches are applied to the decomposition of a set of rolling sub-ranges of yields from the US treasury yields data from the Board of Governors of the Federal Reserve System.

Results

Eigenvector sign change analysis

The choice of covariance estimator had a small but consistent impact on the reliability of PCA in generating eigenvectors that are consistent with the sign change assumptions. This can be seen in Table 1, where the constant correlation shrinkage covariance estimator consistently outperformed the sample covariance estimator while the simple shrinkage estimator consistently underperformed.

The Ledoit-Wolf estimator performed identically to the sample covariance estimator, suggesting that the optimiser had (correctly) inferred that there was little to no benefit in shrinkage towards the diagonal matrix, as in the simple shrinkage estimator.

The other, more intuitive pattern in Table 1 is that performance across all estimators improved with longer sample windows, illustrating how the decomposition is unreliable over shorter time frames. The patterns in the performance of the estimators can also be explored over time, by plotting whether or not each estimation is valid on a rolling basis through the dataset, as in Figure 3 and Figure 4. These plots how how the errors in the estimation are highly correlated across estimators and across sample size ranges, with most of the estimators performing poorly on the same time periods.

These results raise a question: Is the poor performance on small sample sizes due to the short time period sampled, or the small number of samples? This can be investigated by re-running the analysis with the same sample sizes, but with quarterly rather than monthly data, so each slice contains the same number of observations but three times the elapsed time. Comparing these quarterly results in Table 2 to the monthly results in Table 1 suggests that the poor performance on short time windows is largely due to the length of time rather than the number of samples, likely due to insufficient variation in yields over the shorter periods of one-to-two years to provide the required information.

Visual Inspection

We can use visual inspection to understand the results of the sign change analysis in more detail and to find other inconsistencies in the PCA results that were not found by the sign change measure.

Figure 5 and Figure 6 show examples of eigenvectors in one of the sample periods that the constant correlation shrinkage estimator produced valid sign changes but the sample covariance matrix did not. We can see that the cause of the issue is the 3 sign changes in the curvature vector in the sample covariance graph.

We can also see some interesting ways that the eigenvectors can be correct from a sign change perspective, but lack the properties to make them useful in practice. The sample covariance matrix tended to produce very flat eigenvectors when estimated on short time spans of data, as in Figure 7. This lack of structural difference between the eigenvectors after the first few months of maturity limits their usefulness, as they suggest near equal impacts from the principle components on medium and long term bonds.

Conclusion

The main findings are:

Decompositions are unreliable over short periods: Short time periods do not provide the required information to estimate the principle component loadings reliably.

Time period matters more than samples: Poor results over short time periods seem attributable to a lack of time, rather than the small sample sizes. The same small sample sizes spread over
a longer period performed much better, and practically as well as more frequent samples over the same longer period.

The constant correlation shrinkage estimator performs well: The constant correlation shrink- age estimator consistently outperformed the other estimators, including the sample covariance matrix.

The simple shrinkage estimator performs poorly: The simple shrinkage estimator consistently underperformed the sample covariance matrix. This lead to other estimators such as the Ledoit-Wolf estimator, which also shrinks between the sample and diagonal matrix, to pre- form poorly as well, with limited shrinkage benefits.

Covariance matrices are pivotal in finance, serving as a key input for portfolio management, risk assessment, and asset pricing. In portfolio construction, understanding the relationships between different assets is essential for achieving diversification and minimising overall risk. Modern Portfolio Theory relies on covariance matrices to optimize portfolios based on the trade-off between expected return and volatility. Covariance matrices are also used in asset pricing models, aid in pricing derivatives accurately, and for risk management, facilitating the measurement of portfolio risk through metrics like VaR and CVaR. Given the wide range of models that require covariance matrices as an input and the sensitivity of the models to the input matrices, their estimation methods deserve some interrogation.

The obvious (or naive) choice to estimate covariance is the sample covariance matrix. This has a number of key issues that can make it unsuitable for use, including:

Error: While the estimator is unbiased, it is based on a finite sample and so does exhibit estimation error. This error can be significant, particularly when dealing with limited sample sizes, and it affects the accuracy of statistical inferences and predictions based on the estimated covariance matrix.

Sensitivity to noise: The sample covariance matrix is highly sensitive to outliers or extreme values in the data. Outliers can have a disproportionate impact on the estimated covariance, leading to distorted results

Singularity and Inversion: When the number of variables is equal to or greater than the number of observations, the sample covariance matrix becomes singular (non-invertible).

This article will deal primarily with the first two issues of error and sensitivity.

Overview of available estimators

There are a range alternative approaches to covariance matrix estimation

Simplified sample estimators: These methods aim to produce a better matrix by applying sim- plifying assumptions. These include the diagonal covariance matrix, which assumes all vari- ables are independent, and the constant correlation matrix, which assumes all variables have the same correlation.

Shrinkage estimators: Shrinkage estimators aim to produce more useful covariance matrices by combining 2 estimates: A unbiased but high error estimate, F , usually the sample covariance matrix, and a biased but low error estimate, S. These are combined in a weighted sum.

The shrinkage amount δ can be either a given input or be optimised by the estimator.

Time series estimators: As finacial data are often time series, these models aim to incorporate the information provided by the relationships between time periods into the estimated covari- ance matrix. The methods include the Exponentially Weighted Moving Averages (EWMA) model and GARCH (Generalized Autoregressive Conditional Heteroskedasticity) model.

Sparse covariance estimators: When dealing with high-dimensional data, methods like graph- ical lasso (L1 regularization) are used to estimate a sparse covariance matrix, assuming that most of the correlations are zero.

Maximum likelihood estimators: The Maximum Likelihood Covariance Estimator (MLCE) is an approach that seeks to find the covariance matrix that maximizes the log likelihood of the observed data. As with the sample covariance matrix, MLCE is an unbiased estimator, but is known to be sensitive to outliers in small sample sizes and estimation error when the number of variables exceeds the number of observations.

This article will focus on simplified sample estimators and shrinkage estimators.

Quantitative comparison

The impact of estimator choice is evaluated in the context of the use of the estimated matrix in optimising a portfolio of individual equities. This measured through a highly simplified simulation of an investing process with the following steps:

1. The investor accesses return data for their given asset universe, until the current as-at date.
2. The investor calculates the expected returns and and the covariance matrix using the method being tested.
3. The investor determines their portfolio weights using a simple optimiser, with the objective of maximising Sharpe ratio.
4. The as-at date is advanced by the simulation time step of one month, with the portfolio value incriminating by the realised return for the period.

This simulation is run for each of the in-scope covariance measures and for a range of asset universes, being the top-n members of the SP500. Using the current largest members of an index as a investing universe naturally biases the analysis, but this bias is assumed to affect all of the covariance measures equally and so can be ignored.

Methods to be analysed

Sample covariance matrix

The sample covariance matrix is also tested and will act as the baseline that the other methods can be compared against. This is given by the equation below, where X is a matrix of row oriented observations and X is the sample mean vector across all observations.

This is implemented by the python code below

def sample_covariance_matrix(X):
    n = X.shape[1] - 1 # degrees of freedom adj
    avg = np.average(X, axis=1) # Get average of each column
    X = X - avg[:,np.newaxis] # Center data on 0.
    c = np.dot(X, X.T) # Get sum of elementwise differences
    c *= np.true_divide(1, n)
    return c

Diagonal Covariance Matrix

The diagonal covariance matrix simply sets all the non-diagonal elements of the sample covariance matrix to zero. This has the effect of assuming that all variables are independent, with insignificant correlations. Variance (diagonal) elements are left unchanged.

def diagonal_covariance_matrix(X):
    F = sample_covariance_matrix(X)
    n = F.shape[0]
    return F * np.identity(n)

Constant Correlation Covariance Matrix

The constant correlation covariance matrix re-scales each non-diagonal element such that every pair of variables have the same correlation. The chosen correlation is the average correlation, calculated as the average of the non-diagonal elements of the sample correlation matrix. Variance (diagonal) elements are left unchanged.

def constant_correlation_estimator(X):
    F = sample_covariance_matrix(X)
    n = F.shape[0]
    # Calculate correlation
    corr_ = np.ones((n,n))
    for i in range(n):
        for j in range(n):
            corr_[i,j] = F[i,j] / (np.sqrt(F[i,i] * F[j,j]))
    # Get non-diagonal elements
    indicies = np.where(~np.eye(n,dtype=bool))
    avg_corr = np.mean(corr_[indicies])
    # Scale covariances
    S = np.ones((n,n))
    for i in range(n):
        for j in range(n):
            if i != j: # Scale off-diagonal componants
                k = avg_corr / corr_[i,j]
                S[i,j] = F[i,j] * k
            else:
                S[i,j] = F[i,j] # Leave diagonal unchanged
return S

Simple Shrinkage Estimator

The simple shrinkage estimator interpolates linearly between the values of the values of the sample covariance matrix F and the diagonal covariance matrix S.

A shrinkage factor δ of 0.5 is used.

def simple_shrinkage_estimator(X, shrinkage = 0.5):
    F = sample_covariance_matrix(X)
    S = diagonal_covariance_matrix(X)
    return (1-shrinkage)*F + shrinkage*S

Constant correlation shrinkage estimator

This estimator is similar to the simple shrinkage estimator, but rather than shrinking towards the diagonal covariance matrix, it uses the constant correlation matrix.

Ledoit-Wolf Shrinkage Estimator

The Ledoit-Wolf estimator is an extension of the simple shrinkage estimator where an optimal shrinkage factor is estimated based on the input data. The estimator is designed to be computationally efficient and provides a robust estimate, especially when the number of variables is large relative to the number of observations.

The scikit-learn implementation is used for simplicity.

from sklearn.covariance import LedoitWolf
def ledoit_wolf_shrinkage_estimator(X):
    cov = LedoitWolf().fit(X.T)
    return cov.covariance_

Oracle Shrinkage Estimator

Similar to the Ledoit-Wolf estimator, the Oracle estimator extends the simple shrinkage estimator by estimating an optimal shrinkage factor. While the Ledoit-Wolf method typically uses a diagonal matrix as the target, assuming most variables are uncorrelated, the Oracle method allows for more flexibility by incorporating a target matrix derived from prior information or assumptions about the data structure. The Oracle Covariance Estimator also explicitly optimizes the shrinkage parameter to minimize the mean squared error, taking into account the target matrix. These extensions make it more computationally intensive that the Ledoit-Wolf method.

from sklearn.covariance import OAS
def ledoit_wolf_shrinkage_estimator(X):
    cov = OAS().fit(X.T)
    return cov.covariance_

Out of scope methods

The below methods were not analysed, but would be interesting to include in any further analysis.

Factor Models:

Factor models are a class of models used in covariance matrix estimation that assume the observed variables can be decomposed into a linear combination of common factors and idiosyncratic, or specific, factors. The factor model is described as a linear regression of the assets X on their means μ and the common factors F . The residuals U are the idiosyncratic factors

The covariance matrix of the observed variables X can be expressed as:

• Σ is the covariance matrix of X.
• Λ is a diagonal matrix containing the variances of the common factors.
• Ψ is a diagonal matrix containing the variances of the idiosyncratic factors.

If the factors are assumed to completely explain the relationships between the assets Ψ and U can be assumed to be zero.

Results

The impact of the covariance measures are measured through their impact of portfolio return. The time series of portfolio returns for the N asset universe show that the impact on return from the covariance measures are small relative to the overall market movements, as in Figure 1.
Analysing aggregate returns for all covariance measures and asset universe in Table 1 shows that despite returns being highly correlated between the covariance measures, there are significant dif- ferences between aggregate returns for a given asset universe.

Comparing aggregate returns to the sample covariance baseline provides the following insights:

• The Ledoit-Wolf shrinkage estimator provides superior performance over almost all asset universe sizes.
• The simplified diagonal and constant correlation estimators do not provide increased performance.
• When combined with the sample covariance matrix through shrinkage, the diagonal and constant correlation estimators provide better performance when the number of assets is high.
• All shrinkage estimators provide better performance on average when the number of assets is larger. This could be due to the error in the sample covariance matrix increasing with the number of assets, leading to larger gains through the error reduction that shrinkage provides

The relative performance of each covariance estimator compared to the sample covariance matrix can be seen in the time series of the difference in portfolio values between each measure and the sample covariance matrix in Figure 3 and Figure 2. The time series show how each covariance estimator performed well or poorly with fairly high consistency over the simulation.

The consistency of over or under-performance of each covariance estimator can be explored further by looking at the proportion of simulation time steps or investing decisions where each measure outperformed the sample covariance matrix, as in Table 2. This analysis shows that the Ledoit- Wolf and Oracle estimators outperformed the baseline in the majority of the decisions, over all asset universe sizes. The simple shrinkage estimator and the constant correlation estimator were less consistent and as in the aggregate returns analysis performed better on larger asset universe sizes.

Conclusions

There are a number of simple conclusions that can be drawn from these limited results, in the context of the usage of the covariance matrix in portfolio optimisation:

Sample covariance is usually inferior: The sample covariance matrix generally provided worse performance than the shrinkage estimators.

Shrinkage targets alone are not effective: The shrinkage targets themselves, the diagonal and constant correlation matrices, did not provide reliable benefits over the sample covariance matrix

Ledoit-Wolf and oracle estimators are effective: Ledoit-Wolf and oracle estimators provided fairly consistent out performance compared to the sample covariance matrix.

Areas for further research

Factor Models: Comparing the factor model covariance matrix to the estimators explored here.

Multi-asset cases: Asset allocation over many asset classes provides a similar portfolio optimi- sation use case, but with an increase importance on estimating negative correlations, for example between equities and bonds. This will likely lead to measures such as the constant correlation estimator to significantly under-perform.

Other use cases: Covariance matrices are used in many other applications outside of portfolio optimisation, which provide other avenues to investigate the impact of using different estima- tors.

Subscribe to DDIntel Here.

Have a unique story to share? Submit to DDIntel here.

Join our creator ecosystem here.

DDIntel captures the more notable pieces from our main site and our popular DDI Medium publication. Check us out for more insightful work from our community.

DDI Official Telegram Channel: https://t.me/+tafUp6ecEys4YjQ1

Follow us on LinkedIn, Twitter, YouTube, and Facebook.

Perhaps don’t keep-it-super-simple was originally published in DataDrivenInvestor on Medium, where people are continuing the conversation by highlighting and responding to this story.

Motivation

Spotify is a fantastic service which has drastically changed how we find, listen and share music. Their tech forward approach has also allowed new features to be released that further elevate their user’s ability to interact with their music, such as the discover weekly playlist or the recently released blend tool that merges two people’s music tastes into a single playlist.

At the heart of the Spotify user experience is the playlist, which allows the user a way to catalogue and organise songs in a manner reminiscent of the cassette mixtapes of the 1980s. Without playlists the listening experience on Spotify would be a transient one, with users needing to search through the endless library of tracks to find “their” music. As such playlists provide the user with some level of ownership and control over what is “their” music in a system with no ownership and a universe of music and genres.

Given the importance of playlists and the mechanism by which new songs are added to a playlist (just slap them on the end), it is unsurprising that personal playlists can quickly become incoherent, a consequence of evolving tastes and the level of effort required to reorder and reorganise songs.

To avoid this fate for my own playlists, I have started to build solutions to this issue through using Spotify’s developer tools to create my own methods to interact with the Spotify music universe.

This is part one in a series, with links to the other parts at the bottom of this article.

The Spotify web API

Spotify offers a free Web API that allows you to programmatically access various aspects of Spotify and modify parts of a user’s account — given their permission at least!

I want to be able to access a user’s playlists and tracks within, as well as be able to modify these playlists or create new ones. This requires 4 key abilities:

• Ability to list all of a users playlists
• Ability to list the tracks within a playlist
• Ability to get data about a track
• Ability to create or modify a user’s playlist

Fortunately the web API offers all of these functions (and many more), as well as the ability to bundle many requests into a single HTTPS call, which reduces the number of requests required when loading track data on all tracks within a playlist.

For reference, the required playlist level interactions are contained within the Playlist API while the track level interactions are contained within the Tracks API.

The basic flow can be summarised as below:

Authorisation

For any of the API features to be useful, the user must be logged in as a valid Spotify user. Spotify provides a number of ways for a user to log in and provide the application with the required authorisation to access their data and make changes.

For my applications I have chosen to use the implicit grant method as I do not need long-term persistent authorisation and as implicit grant authorisation allows the application to be front-end only. This flow provides the application with a token that is then used in all subsequent API calls.

Pagination

Due to the large number of playlists and songs a user may have, all of the used API calls have a maximum return length as well as pagination to get the next chunk of data if any exists. This keeps the size of each api call manageable, especially when the data is contained within the URL parameters rather than headers, and so is subject to the strict character limits by browser.

The tool

For those that are interested the playlist optimiser tool is available at https://www.grantholtes.com/smartshuffle

Other stories

• Part 2: Playlist Optimisation as a Traveling Salesperson Problem

Building Better Playlists: Spotify Integration was originally published in CodeX on Medium, where people are continuing the conversation by highlighting and responding to this story.

Motivation

Spotify is a fantastic service which has drastically changed how we find, listen and share music. Their tech forward approach has also allowed new features to be released that further elevate their user’s ability to interact with their music, such as the discover weekly playlist or the recently released blend tool that merges two people’s music tastes into a single playlist.

At the heart of the Spotify user experience is the playlist, which allows the user a way to catalogue and organise songs in a manner reminiscent of the cassette mixtapes of the 1980s. Without playlists the listening experience on Spotify would be a transient one, with users needing to search through the endless library of tracks to find “their” music. As such playlists provide the user with some level of ownership and control over what is “their” music in a system with no ownership and a universe of music and genres.

Given the importance of playlists and the mechanism by which new songs are added to a playlist (just slap them on the end), it is unsurprising that personal playlists can quickly become incoherent, a consequence of evolving tastes and the level of effort required to reorder and reorganise songs.

To avoid this fate for my own playlists, I have started to build solutions to this issue through using Spotify’s developer tools to create my own methods to interact with the Spotify music universe.

This is part two in a series, with links to the other parts at the bottom of this article.

Playlist Optimisation

“DJing is a complex analytical as well as creative task, creating an automatic DJ has proved to be highly nontrivial” Len Vande Veire and Tijl De Bie, EURASIP Journal on Audio, Speech, and Music Processing (2018)

User’s Spotify playlists tend to be chronologically organised by the date when a song is added, which does not typically lead to a well ordered playback. When seeking to improve this order, it helps to reframe the problem as “what order should we play the songs such that every subsequent song is acceptably similar to the previous song”.

One shortcut to determine how similar two songs is to use the Spotify audio analysis data, which includes measurements of tempo, key, as well as more subjective measures such as “danceability” and “acousticness” which capture the sound and mood of the song. If a set of these measures are stored as vectors for each song, the similarity between songs can be quantified as the distance between these vectors, with various dimensions weighted by relative importance. Two songs that are nearby as measure by their vectors should have a similar BPM, key and general mood, while songs that are further apart may be very different and so unsuitable to play back-to-back.

The task of ordering the songs of the playlist to minimise the sum of the distance between sequential songs is version of the traveling salesman problem or TSP, which seeks find a route between cities that minimise the total distance traveled. The only major difference is that instead of a two dimensional real-world space, we have an arbitrary number of dimensions in the song vectors.

Finding an exact solution to the TSP is a famously NP-hard problem, which requires an infeasible amount of computational power for larger numbers of data points. The complexity of the TSP is O(n!), which becomes impractical beyond roughly 20 data points, and even optimised versions are impractical beyond 50 data points. As we don’t know how long a user’s playlist may be and want the algorithm to run in a standard web browser, seeking an exact solution is unacceptable.

Instead an approximate solution is used, in this case a simple nearest neighbours approach is used which has a complexity of O(n squared). The algorithm can be summarised as the following steps.

1. Mark all points as unseen.
2. Select a point from the set, in this case the array of playlist vectors. Append the point to the output array and mark it as seen.
3. While there are unseen points, find the nearest point by looping over all unseen points and measuring the distance to the previously seen point.
4. Append the nearest point to the output array and mark it as seen.
5. Return the output array

This greedy approach is far from optimal, but is significantly faster and simpler than competing approaches, and is well suited to this application where speed is preferred over an exact solution. An additional perk of this method is that every starting point yields a slightly different approximate solution, which allows the user to shuffle between approximately optimal solutions by simply recomputing the solution, given that the starting point is randomly chosen. When this is combined with a visualisation of the solution, it allows the user to pick the solution they like the look of best. This solution could be further improved by selecting a starting point with a desirable characteristics for the start of a playlist, such as low BPM compared to the other songs.

The tool

For those that are interested the playlist optimiser tool is available at https://www.grantholtes.com/smartshuffle

Other stories

• Part 1: Building Better Playlists: Spotify Integration

Playlist Optimisation as a Traveling Salesperson Problem was originally published in DataSeries on Medium, where people are continuing the conversation by highlighting and responding to this story.

For a recent project I needed to find a method to determine which colour palettes were present in a large dataset of artworks, ideally without any input from myself! This task aligned itself to the unsupervised technique of clustering and can be divided into two steps:

1. Extracting commonly occurring colours.
2. Extracting common combinations of these colours.

Both of these can be accomplished with K-Means clustering, which simply clusters points by the distance between points and cluster centers.

1. Extracting commonly used colours.

To effectively use K-means to extract colour combinations it is important to have the euclidean distance between colour data points be interpretable as perceived difference in colours. Simply using the red-green-blue (RGB) channels of the images as a three dimensional vector does a rather poor job of this, as distance in RGB space does not translate well into perceived difference.

To overcome this the image is converted into YCbCr colour space, with the Y, Cb, Cr channels forming three-dimensional vectors to represent each color. Distances in this colour space do a much better job at representing differences in perceived colour.

After converting the loaded dataset into YCbCr and extracting a large random sample of pixels, K-means with n clusters will extract n commonly used colours as the cluster centers. As these centres mark the average of the colour clusters, the colours produced can be rather muted if an insufficient number of clusters are used. Empirically I found that 20 to 40 clusters produced good results.

2. Extracting common combinations of these colours.

Once the common colours have been found, we want to find out in what combination and proportions these colours are commonly used.

This can be done by taking a fixed sample of pixels from each image and seeing which colour cluster each pixel belongs to. These groupings can be aggregated as a count of each colour by image and stored as an n-dimensional vector that represents the colour usage in each image. Using counts in this manner captures the proportion as well as presence of each colour. As all of these image vectors have the same sum, these image representations lie on the hyperplane that intersects all axes at the size of the pixel sample.

K-means can then be applied to these count vectors to extract the common combinations and proportions of the colours found before, with m-cluster centres yielding m different colour palettes.

This system can be tested by providing a dataset that represents images from two very different colour pallete distributions — bright rainbows and earthy desserts. Setting the number of colours at a conservative 15 and the number of colour pallets to two we can see how the different palettes are extracted through the two steps.

Extracting Colour Palettes with Unsupervised Learning was originally published in Analytics Vidhya on Medium, where people are continuing the conversation by highlighting and responding to this story.

Pandemics, bushfires and economic crisis. Not only the menu for the last year but a small subset of the phenomena that can be modelled under the framework of agent based modelling (ABM).

ABM is well suited to understanding how complex behavior emerges in systems based on simple interactions between the system’s individual participants. ABMs have shown particular strength in explaining general observations, referred to as ‘stylized facts’, such as observed neutron scattering patterns or the distribution of asset returns. These agents could be anything from people or companies to animals in an ecosystem or atoms in a gas. This approach to modelling is a bottom up approach, where the system’s behavior is completely driven by assumptions on participants’ individual behavior, rather than a top down approach which assumes on how the participants will behave together.

A simple example of this is wildfire modelling, with each agent representing a small area of land. The rules for each patch of land can be extremely simple: If the land is on fire, set downwind neighbouring patches of land alight until all fuel in the patch is burnt. Despite these simple rules, the behavior of the wildfire system as a whole displays complex behavior. Furthermore, exploring the effects of changing assumptions such as making fire spreading probabilistic or introducing fuel regrowth only requires the rules around interactions between land patches to be altered, leaving the model framework unchanged. This aids in rapid experimentation and testing of assumptions.

Rapid experimentation isn’t the only advantage of ABM over top down approaches. Accounting for variation in the individuals being modelled can introduce significant complexity in a top down model while being a trivial change in an agent based approach, as this heterogeneity can be captured by varying individual agents. This can be used to capture differences such as varying fuel levels in wildfire simulations or create populations of agents that mimic the census data for epidemic models.

The ABMs are not limited to agents interacting in simulated physical space, with the model framework extending to non-physical limitations on agent interaction such as limiting agents to only interacting with others that they ‘know’. This allows for social heterogeneity to be introduced, with many agents belonging to closed social circles while fewer highly social agents bridge these circles.

This is achieved in the simple epidemic model below through expressing the agents social connections as a network with small-world properties.

This also allows us to explore the stylized facts around how limiting social interaction affects the duration and intensity of the epidemic. In the model this is achieved by reducing the average connections in the network, which is akin to individuals reducing the physical contacts to only their closest friends and family.

Another key advantage of ABM is the explainability of the underlying assumptions, which can provide confidence in the results and make findings easier to communicate. For example, the top down analysis of income inequality and the impact of taxation is usually highly mathematical and very abstracted from individuals, making communication of policy decisions a consistent challenge.

On the other hand, a bottom up ABM approach allows for the model to be reduced to only three steps, which are repeated by agents each cycle: 1) randomly trade with agents that they have social links to, as in the epidemic model; 2) pay tax on income from trades and receive welfare payments if wealth is below a minimum level; 3) incur a fixed living expense and receive interest on savings. Noting that these assumptions say nothing about the work ethic or productivity of each agent, the stylized facts from the model still emulate those found in top down models, particularly that flat tax results in much higher wealth inequality (measured as the Gini coefficient) than progressive tax, and that welfare payments are required to prevent fixed living costs from bankrupting a large proportion of the agents.

Unfortunately there is often a lack of knowledge around these underlying assumptions and inputs, as relevant data on the systems in question is usually captured at an aggregate level rather than an individual level. Furthermore, the model outputs are often sensitive to the inputs used which can greatly limit use in forecasting and erode confidence.

One solution is to capture the assumptions as a set of parameters, then optimise the parameters to fit historical data, running the model for various parameter combinations until the model outputs match the observed data. The parameters from this tuned model are then used for ongoing forecasting or modelling exercises. This practice increases confidence in the insights generated and provides stronger links between the simulated and real worlds.

While ABM techniques do not replace traditional top down approaches, they do offer an alternative perspective and in some cases provide a method to approach problems that would be too complex to define in a top down manner or where the agents display a great deal of heterogeneity. While there is a level of technical skill required to write and test ABMs, this learning curve flattens due to the reusability of the model framework and the existence of open source libraries such as Mesa for Python which provide the general framework for agents, agent interactions and data capture. Finally, confidence in model results is increased through simple-to-communicate assumptions and methodologies. This is further tested and refined by basing model inputs on known real world data and optimising parameters to fit historical observations.

This story is the first in a series on agent based modelling, and is based on my own research and experimentation with building ABMs. If you want to follow along please consider following my profile or this publication.

Some sources are below:

https://www.sciencedirect.com/science/article/pii/S1755436517300221

https://www.bankofengland.co.uk/-/media/boe/files/quarterly-bulletin/2016/agent-based-models-understanding-the-economy-from-the-bottom-up.pdf

A single technique to model forest fires, pandemics and economies was originally published in TDS Archive on Medium, where people are continuing the conversation by highlighting and responding to this story.

Colour combination analysis with chained clustering

For a recent project I needed to find a method to determine which colour palettes were present in a large dataset of artworks, ideally without any input from myself! This task aligned itself to the unsupervised technique of clustering and can be divided into two steps:

1. Extracting commonly occurring colours.
2. Extracting common combinations of these colours.

Both of these can be accomplished with K-Means clustering, which simply clusters points by the distance between points and cluster centers.

1. Extracting commonly used colours.

To effectively use K-means to extract colour combinations it is important to have the euclidean distance between colour data points be interpretable as perceived difference in colours. Simply using the red-green-blue (RGB) channels of the images as a three dimensional vector does a rather poor job of this, as distance in RGB space does not translate well into perceived difference.

To overcome this the image is converted into YCbCr colour space, with the Y, Cb, Cr channels forming three-dimensional vectors to represent each color. Distances in this colour space do a much better job at representing differences in perceived colour.

After converting the loaded dataset into YCbCr and extracting a large random sample of pixels, K-means with n clusters will extract n commonly used colours as the cluster centers. As these centres mark the average of the colour clusters, the colours produced can be rather muted if an insufficient number of clusters are used. Empirically I found that 20 to 40 clusters produced good results.

2. Extracting common combinations of these colours.

Once the common colours have been found, we want to find out in what combination and proportions these colours are commonly used.

This can be done by taking a fixed sample of pixels from each image and seeing which colour cluster each pixel belongs to. These groupings can be aggregated as a count of each colour by image and stored as an n-dimensional vector that represents the colour usage in each image. Using counts in this manner captures the proportion as well as presence of each colour. As all of these image vectors have the same sum, these image representations lie on the hyperplane that intersects all axes at the size of the pixel sample.

K-means can then be applied to these count vectors to extract the common combinations and proportions of the colours found before, with m-cluster centres yielding m different colour palettes.

This system can be tested by providing a dataset that represents images from two very different colour pallete distributions — bright rainbows and earthy desserts. Setting the number of colours at a conservative 15 and the number of colour pallets to two we can see how the different palettes are extracted through the two steps.

Extracting Colour Palettes with Unsupervised Learning was originally published in DataSeries on Medium, where people are continuing the conversation by highlighting and responding to this story.

The range and accessibility of basic services such as supermarkets and banks varies along with average income between neighboring postcodes in Victoria, Australia.

However, unlike the US where “banking deserts” are leaving low-income communities without access to financial services and western Sydney where “food deserts” are limiting access to groceries, Victoria seems to demonstrate the opposite pattern.

Lower income communities have significantly more schools, banks, supermarkets and health vendors within walking distance than their high income counterparts. Conversely, high income suburbs do have significantly more options within a larger driving distance.

• A high income suburb has 40% fewer banking institutions and 51% fewer grocery stores within walking distance than neighboring low income suburbs.
• Higher income suburbs have 48% more health providers and 30% more schools within driving distance.

This somewhat counter intuitive result suggests that basic services may be more numerous and distributed more widely in and around higher income suburbs, while lower income suburbs enjoy fewer, more densely grouped and accessible services.

These results are based on an analysis of the ATO postcode income data and location data extracted from the Google Maps API.

JsonMapsData = json.loads( requests.get("https://maps.googleapis.com/maps/api/place/textsearch/json?query=supermarkets+in+3052&radius=10000&key=YOUR_API_KEY).text)

To account for differences between regions within the state, postcodes are compared to other nearby postcodes within subregions. To further control for intra-subregion differences in general accessibility, the distance to the nearest regional center and urban or regional status of the postcode are also used as controls. Poisson regression was used to provide estimates of the intra-subregion variation in service availability. For the technically inclined, the code and data is available here.

The Surprising Pattern of Food and Service Deserts in Victoria was originally published in DataSeries on Medium, where people are continuing the conversation by highlighting and responding to this story.

In the previous story, I developed a model to encode the art style of an image into a vector in higher dimensional space, where the euclidean distance between vectors represents how visually similar the images are.

For an image recommendation system to use these embeddings, the vector for every image in the data set needs to be stored on disk and the distance between two or more vectors will need to calculated at some point. The more dimensions that are present in each vector, the greater the storage and computational needs of the system, increasing costs and search times.

Due to this, it is beneficial to reduce the number of dimensions in the embedding vector, while minimizing the information lost in doing so. To achieve this a number of techniques were tested, and the methods used are explored in more detail below. However, for these embeddings the more traditional methods outperformed neural network based methods. Due to this, PCA dimensionality reduction is used moving forward, with the PCA trained on a sample of the embeddings before being applied to the remainder in order to reduce the time taken to train the PCA.

Efficiency of dimensionality reduction

Efficient dimensionality reduction from the perspective of information relies on at least some of the dimensions present being redundant, meaning that they can be replaced by a function of the other dimensions present.

For example, consider a 3D vector V = [X, Y, Z]. If Z can be expressed as an exact function of Y, X it can be removed from the vector without removing any information from the system. We can always reconstruct Z later by knowing X and Y. Even if Z is an inexact function of X, Y, for example Z = X + Y + e where e is an unobserved variable, we can still remove Z at the cost of removing the information provided by e. The information lost increases as the variance of e increases. Overall, the information lost through dimensionality reduction can be found by reducing the dimensions of V, then restoring the removed dimensions and comparing the variance of the original and transformed vectors. The lost variance is a metric for the information lost through dimensionality reduction.

The embedding vectors I am working with have characteristics which make them good candidates for dimensionality reduction, consisting of two parts:

1. A (flattened) co-variance matrix between textures: This matrix provided the co-variance between every combination of two textures. Consider three textures i, j, k. If textures i and j are correlated and j and k are correlated, we would expect i and k to be correlated. This means that the dimension for the co-variance between i and k contains some redundant information. The convenience matrix contains thousands of dimensions, many of which will be redundant and can be removed without sacrificing a large amount of information.
2. A vector of what proportion of the image is made up of each colour, which each dimension representing a colour. Some colours may never be used so can be removed, while other colours may show high enough correlations to be redundant.

PCA

PCA is a linear dimensionality reduction technique, exploiting linear relationships between variables to express the information they hold in fewer dimensions. This is shown in the example below, where two dimensions are condensed into one.

This technique was very efficient, reducing the dimensions by 94% while retaining 99.3% of the variance of the original data.

Neural Network Autoencoders

A variety of neural network based autoencoders were also tested. These take the embedding vectors as both the input and the target, and use successive layers to between to decrease the size of the embedding down to a small “encoded” state then back up to its original dimensions in an hourglass like fashion, as below.

Much like PCA, this setup minimises the error between the input and the decoded vectors, while creating an intermediate encoded vector with far fewer dimensions. Unlike PCA, this technique introduces non-linearities through multiple layers, allowing for more complex relationships between the elements of the embedding vector to be captured in the lower dimension representation. Ideally, this should result in a lower amount of information loss, however, this technique resulted in significantly more information loss and longer compute times than PCA.

Superior information loss may be possible through further experimentation with different network architectures, but the accuracy and quick compute times offered by PCA made it an obvious choice out of the techniques tested.

Comparison of Dimensionality Reduction Techniques for Image Art Style Embeddings was originally published in TDS Archive on Medium, where people are continuing the conversation by highlighting and responding to this story.

The aim of this exercise is to find a function that transforms images into embedding vectors where the euclidean distance between vectors represents how visually similar the images are. This allows a nearest neighbors search on one image’s embedding to return images that are visually similar, empowering image recommendation and clustering.

This task is complicated by the images being illustrations rather than photographs. While photographs are true representation of the subject, and hence similar subjects and composition would yield visually similar images, illustrations include the extra dimension of style.

Style: “a way of painting, writing, composing, building, etc., characteristic of a particular period, place, person, or movement.”

To create the image embeddings two methods are used:

• Convolutional Neural Network embeddings
• 3D colourspace nearest neighbours

Convolutional Neural Network embeddings

The model used to transfer images to embeddings borrows findings from both neural-style transfer and ecommerce image content similarity work. A Neural Algorithm of Artistic Style by Gatys et. al and writing by @Raghul Asokan was used to research neural-style transfer. Creating an image similarity function with TensorFlow and its application in e-commerce by Nina Pakhomova was used to develop an understanding of similar image retrieval techniques.

Both style transfer and image similarity techniques leverage pre-trained convolutional (CNN) image classification models to extract and abstract information from images. These models have been trained to extract information for classification purposes, with the intermediate convolutional layers capturing increasingly complex information about the contents of the image as the depth of the network increases.

For this I will be using VGG-19 model trained on the image-net database, the architecture of which is shown in the diagram above. To extract artistic style using this network we must first make a few subjective assumptions about what defines style in the context of the task at hand.

• We assume that a style can be described as a combination of various textures and colours used together.
• We also assume that sub-features, small repeating parts of the overall image, are important to the style we experience. For example, a pencil drawing of a building and of a landscape my have similar textures, but the sharp right angles and windows of the building lend the architectural drawing a style that is distinct to that of the landscape.

Luckily we know that CNNs capture and extract this information in each convolutional layer, with each kernel “searching” for a texture or sub-feature. When that texture or sub-feature is present in the underlying image, the kernel outputs a high value when it samples that region of the image. The activation process for a single kernel is shown below.

When all kernels for a layer are passed over the input, a 3D matrix with a depth equal to the number of kernels, k is output.

The deeper into the network these kernels are, the more complex the extracted textures, with the deepest kernels activating with entire structures such as eyes or windows. This behavior is demonstrated in deep dream images.

Using our above assumptions about what defines art style and knowing that each kernel captures different distinct elements of art style, we can use the combination of activations from all kernels at a specific region as a mathematical description of the art style that region. To summarize these findings over all regions, one could simply take every single activation and flatten them into an embedding vector, as is done in ecommerce image content similarity. However, that approach is inefficient as it requires that every activation is stored and is sensitive to spacial variations over the image. That would result in the two images below being rated as very dissimilar in style despite only varying spatially.

To overcome these issues Gatys et. al use gram matrices in their neural style transfer paper. Gram matrices compute the correlation between the outputs of every kernel. Kernel combinations that show high correlation in their outputs are taken to be descriptive of the style of the image. For example in Starry Night above, the textures of individual brushstrokes and swirls would be highly correlated over both of the images, and as such the gram matrices of the images would be similar.

Computationally, gram matrices are calculated by flattening the convolutional output into a 2D array, with columns for each kernel and rows for every activation, A. The dot product of this array with its transpose is taken and output as the gram matrix, G. Here G(i,j) gives the correlation between the activations of kernels i and j. Note that the term “correlation” is used rather loosely as this is not the statistically true correlation. Potentially the term cumulative co-activation would be more accurate.

The size of the gram matrix produced is k² which is flattened into a vector for storage. Which intermediate layer(s) to extract the kernel outputs from was determined by trial and error, with emphasis placed on performance, as the deeper the layer(s) used the more calculations and compute time required.

3D colourspace nearest neighbours

Which colours are present is an obvious visual descriptor and one that we did not want to overlook. While gram matrices will extract colour information, we used a secondary sampling approach to explicitly encode the dominant colours present. Colour dimensionality reduction is well documented and usually to extract the dominant colours in an image I would use a clustering approach such a K-means. However, in this case I want the express the colours of the image as a vector where each element always corresponds to the same colour, allowing comparison between images.

To achieve this the 3D RGB colourspace of an image is divided into n³ mutually exclusive cubes, each with side length 255/n. The number of pixels within each cube are counted, and the counts are stored in a vector of length n³.

Embedding construction

Finally, the flattened gram matrices are concatenated with the colour embeddings to form embeddings of dimensionality N = k² + n³.

Results

A sample set of 10,000 images was used to test the model. As this method is unsupervised and requires no further training this small set is sufficient.

If you are reading this, thanks for making it this far! Feel free to follow along as this project continues.

Developing Art Style Embeddings for Visual Similarity Comparison of Artworks was originally published in TDS Archive on Medium, where people are continuing the conversation by highlighting and responding to this story.