DIET NETWORKS: THIN PARAMETERS FOR FAT GENOMICS
Learning tasks such as those involving genomic data often poses a serious challenge: the number of input features can be orders of magnitude larger than the number of training examples, making it difficult to avoid overfitting, even when using the known regularization techniques. We focus here on tasks in which the input is a description of the genetic variation specific to a patient, the single nucleotide polymorphisms (SNPs), yielding millions of ternary inputs. Improving the ability of deep learning to handle such datasets could have an important impact in medical research, more specifically in precision medicine, where highdimensional data regarding a particular patient is used to make predictions of interest. Even though the amount of data for such tasks is increasing, this mismatch between the number of examples and the number of inputs remains a concern. Naive implementations of classifier neural networks involve a huge number of free parameters in their first layer (number of input features times number of hidden units): each input feature is associated with as many parameters as there are hidden units. We propose a novel neural network parametrization which considerably reduces the number of free parameters. It is based on the idea that we can first learn or provide a distributed representation for each input feature (e.g. for each position in the genome where variations are observed in data), and then learn (with another neural network called the parameter prediction network) how to map a feature’s distributed representation (based on the feature’s identity not its value) to the vector of parameters specific to that feature in the classifier neural network (the weights which link the value of the feature to each of the hidden units). This approach views the problem of producing the parameters associated with each feature as a multi-task learning problem. We show experimentally on a population stratification task of interest to medical studies that the proposed approach can significantly reduce both the number of parameters and the error rate of the classifier.
1 INTRODUCTION
Medical datasets often involve a dire imbalance between the number of training examples and the number of input features, especially when genomic information is used as input to the trained pre-
∗Equal contribution.
dictor. This is problematic in the context where we want to apply deep learning (which typically involves large models) to precision medicine, i.e., making patient-specific predictions using a potentially large set of input features to better characterize the patient. This paper proposes a novel approach, called Diet Networks, to reparametrize neural networks to considerably reduce their number of free parameters when the input is very high-dimensional and orders of magnitude larger than the number of training examples.
Genomics is the study of the genetic code encapsulated as DNA in all living organisms’ cells. Genomes contain the instructions to produce and regulate all the functional components needed to guide the development and adaptation of living organisms. In the last decades, advances in genomic technologies resulted in an explosion of available data, making it more interesting to apply advanced machine learning techniques such as deep learning. Learning tasks involving genomic data and already tackled by deep learning include: using Convolutional Neural Networks (CNNs) to learn the functional activity of DNA sequences (Basset package, Kelley et al. (2016), predicting effects of noncoding DNA (DeepSEA, Zhou & Troyanskaya (2015)), investigating the regulatory role of RNA binding proteins in alternative splicing (Alipanahi et al., 2015), inferring gene expression patterns (Chen et al., 2016; Singh et al., 2016) and population genetic parameters (Sheehan & Song, 2016) among others (see Leung et al. (2016) for a detailed example). Noticeably, most of these techniques are based on sequence data where convolutional or recurrent networks are appropriate. When the full DNA sequence is unavailable, such as when data is acquired through genotyping, other methods need to be used. All this work shows that deep learning can be used to tackle genomic-related tasks, paving the road towards a better understanding of the biological impact of DNA variation.
Applying deep learning to human genetic variation holds the promise of identifying individuals at risk for medical conditions. Modern genotyping technologies usually target millions of simple variants across the genome, called single nucleotide polymorphisms (SNPs). These genetic mutations result from substitutions from one nucleotide to another (eg. A to C), where both versions exist within a population. In modern studies, as many as 5 millions SNPs can be acquired for every participant. These datasets differ from other types of genomic data because they focus on the genetic differences between individuals which represents a space of high dimensionality where sequencecontext information is unavailable. In medical genetics, these variants are tested for their association with a trait of interest, an approach termed genome-wide association study (GWAS). This methodology aims at finding genetic variants implicated in disease susceptibility, etiology and treatment.
An important confounding factor in GWAS is population stratification, which arises because both disease prevalence and genetic profiles vary from one population to the other. Although most GWAS have been restricted to homogeneous populations, dimensionality reduction techniques are generally used to account for population-level genetic differences (Price et al., 2006). Our experiments compare such dimensionality reduction techniques (based on principal components analysis, PCA) to the proposed Diet Network parametrization, as well as with standard deep networks.
Recently, several machine learning methods have been successfully applied to detect population stratification, based on the presence of systematic differences in genetic variation between populations. For instance, Support Vector Machines (SVM) models have been used multiple times to infer recent genetic ancestry of sub-continental populations (Haasl et al. (2013)), and local ancestry in admixed populations (SupportMix, Omberg et al. (2012), 23andMe, Inc.). However SVM methods are very sensitive to the the kernel choice and the parameters. They also tend to overfit the model selection criterion which usually induces a limitation in its predictive power.
In this work, we are interested in predicting the genetic ancestry of an individual from their SNP data using a novel deep learning approach, Diet Networks, which allow us to considerably reduce the number of free parameters. Therefore, we propose to tackle this problem by introducing a multi-task architecture in which the problem of predicting the appropriate parameters for each input feature is considered like a task in itself, and the same parameter prediction network is used for all of the hundreds of thousands of input features. This parameter prediction network learns to predict these feature-specific parameters as a function of a distributed representation of the feature identity, or feature embedding. The feature embedding can be learned as part of end-to-end training or using other datasets or a priori knowledge about the features. What is important is that two features which are similar in some appropriate sense (in terms of their interactions with other features or other variables observed in any dataset) end up having similar embeddings, and thus a similar parameter vector as output of the parameter prediction network. A practical advantage of this approach is
that the parameter prediction network can generalize to new features for which there is no labeled training data (without the target to be predicted by the classifier), so long as it is possible to derive an embedding for that feature (for example using just the unlabeled observations of co-occurences of that feature with other features in human genomes).
An interesting consideration is that from the point of the parameter prediction network, each feature is an example: more features now allow to better train the parameter prediction network. It is like if we were considering not the data matrix itself but its transpose. This is actually how the Diet Network implementation processes the data, by using the transpose of the matrix of input values as the input part of the learning task for the parameter prediction network.
The idea of having two networks interacting with each other and with one producing parameters for the other is well rooted in the machine learning literature (Bengio et al., 1991; Schmidhuber, 1992; Gomez & Schmidhuber, 2005; Stanley et al., 2009; Denil et al., 2013; Andrychowicz et al., 2016). Recent efforts in the same direction include works such as (Bertinetto et al., 2016; Brabandere et al., 2016; Ha et al., 2016) that use a network to predict the parameters of a Convolutional Neural Network (CNN). Brabandere et al. (2016) introduce a dynamic filter module that generates network filters conditioned on an input. Bertinetto et al. (2016) propose to learn the parameters of a deep model in one shot, by training a second network to predict the parameters of the first from a single exemplar. Hypernetworks (Ha et al., 2016) explore the idea of using a small network to predict the parameters of another network, training them in an end-to-end fashion. The small network takes as input the feature embedding from the previous layer and learns the parameters of the current layer.
To the best of our knowledge, deep learning has never been used so far to tackle the problem of ancestry prediction based on SNP data. Compared to other approaches that attempt to learn model parameters using a parameter prediction network, our main goal is to reduce the large number of parameters required by the model, by considering the input features themselves as sub-tasks in a multi-task view of the learning problem, as opposed to constructing a model with even higher capacity, as seen, e.g. in (Ha et al., 2016). Our approach is thus based on building an embedding of these tasks (the features) in order to further reduce the number of parameters.
We evaluate our method on a publicly available dataset for ancestry prediction, the 1000 Genomes dataset1, that best represents the human population diversity. Because population-specific differences in disease and drug response are widespread, identifying an individual’s ancestry heritage based on SNP data is a very important task to help detect biological causation and achieve good predictive performance in precision medicine. Most importantly, ancestry-aware approaches in precision genomics will reduce the hidden risks of genetic testing, by preventing spurious diagnosis and ineffective treatment.
2 METHOD
In this section, we describe the Diet Networks as well as the feature embeddings used by the model.
2.1 MODEL
Our model aims at reducing the number of free parameters that a network trained on fat data would typically have.
Let X ∈ RN×Nd be a matrix of data, with N samples and Nd features, where N Nd (e. g. N being approximately 100 times smaller than Nd). We build a multi-layer perceptron (MLP), which takes X as input, computes a hidden representation and outputs a prediction Ŷ. Optionally, the MLP may generate a reconstruction X̂ of the input data from the hidden representation. Figure 1(a) illustrates this basic network architecture. Let xi be one data sample, i.e. a row in X. The standard formulations to compute its hidden representation hi, output prediction ŷi and reconstruction x̂i are given by
hi = f(xi), ŷi = g(hi), x̂i = r(hi), (1)
where f , g and r are non-linear functions.
1http://www.internationalgenome.org/
The number of parameters of the first hidden layer of the architecture grows linearly with the dimensionality of the input data:
h (1) i = f1(xiWe + be), (2)
where We and be are the layer’s parameters. Using fat data such as the one described in Section 1, leads to a parameter explosion in this layer, hereafter referred to as fat hidden layer. To give the reader an intuition, consider the case of having an input with Nd = 300K, and a hidden layer with N1h = 100, the number of parameters of such a layer would be 30M. The same happens to the number of parameters of the optional reconstruction layer, hereafter referred to as fat reconstruction layer.
In order to mitigate this effect, we introduce an auxiliary network to predict the fat layers’ parameters. The auxiliary network takes as input the transposed data matrix XT, extracts a feature embedding and learns a function of this embedding, to be used as parameters of a fat layer:
(We)j: = φ(ej), (3)
where ej represents the embedding of a feature in XT, φ is a non-linear function and (We)j: is the j-th row of We. This means that each feature is associated with the vector of values it takes in the dataset (e.g. across the patients). Other representations could be used, e.g., derived from other datasets in which those features interact. Figure 1(b) shows a prediction network which is an auxiliary network that predicts the parameters of the fat hidden layer of our basic network. Following the same spirit, Figure 1(c) highlights the interaction between a second prediction network that predicts the fat reconstruction layer parameters and the basic network. The architectures of both auxiliary networks may share the initial feature embedding.
The feature embeddings used in the auxiliary networks allow us to substantially reduce the number of free parameters of the fat layers of the basic architecture. The auxiliary network should predict a matrix of weights of size Nd × N1h from a feature embedding. Consider a feature embedding that would transform each N -dimensional feature into a Nf -dimensional vector, where Nf < N . The auxiliary network would learn a function φ : RNf → RN1h . Thus, the fat hidden layer of our basic architecture would have Nf × N1h free parameters (assuming a single layer MLP in the auxiliary network), instead ofNd×N1h . Following our previous example, whereNd = 300K andN1h = 100, using an auxiliary network with previously-obtained feature embeddings of dimensionality Nf = 500 would reduce the number of free parameters of the basic network by a factor of 600 (from 30M to 50K).
The model is trained end-to-end by minimizing the following objective function
H(Ŷ,Y) + γ||X̂−X||22, (4)
where H refers to the cross-entropy, Y to the true classification labels and γ is a tunable parameter to balance the supervised and the reconstruction losses.
2.2 FEATURE EMBEDDINGS
The feature embeddings used by the auxiliary networks can be either pre-computed or learnt offline, as well as learnt jointly with the rest of the architecture. In theory, any kind of embedding could be used, as long as we keep in mind that the goal is to reduce the number of free parameters of the basic model. In this work, we considered random projections (Bingham & Mannila, 2001), histograms (which are akin to bag-of-words representations), feature embeddings learnt offline (Mikolov et al., 2013) and feature embeddings jointly learnt with the rest of the proposed architecture.
Random projection: Randomly initializing an MLP defines a random projection. By using such a projection to encode the high-dimensional feature space into a more manageable lower-dimensional space, we were able to obtain decent results.
Per class histogram: For a given SNP, we can define a histogram of the values it can take over the whole population. Once normalized, this yields 3 values per SNP, corresponding to the proportion of the population having the values 0, 1 and 2 respectively for that SNP. After initial tests showed this was too coarse a representation for the dataset, we instead chose to consider the per-class proportion of the three values. With 26 classes in the 1000 Genomes dataset, this yields an embedding of size 78 for each feature. By this method, the matrix XT is summarized as a Nd × 78 matrix, where Nd is the number of SNPs in the dataset.
SNPtoVec: In Mikolov et al. (2013), the authors propose a word embedding that allows good reconstruction of the words’ context (surrounding words) by a neural network. SNPs do not have a similarly well-defined positional context (SNPs close together in our ordering might very well be independent) so our embedding is instead built by training a denoising autoencoder (DAE) (Vincent et al., 2008) on the matrix X. Thus, the DAE learns to recover the values of missing SNPs by leveraging their similarities and cooccurences with other SNPs. Once the DAE is trained, we obtain an encoding for each feature by feeding to the DAE an input where only that feature is active (the other features are set to 0s) and computing the hidden representation of the autoencoder for that single-feature input.
Embedding learnt end-to-end from raw data: In this case, we consider the feature embedding to be another MLP, whose input corresponds to the values that a SNP takes for each of the training samples and, whose parameters are learnt jointly with the rest of the network. Note that the layer(s) corresponding to the feature embedding are shared among auxiliary networks. For experiments reported in Section 4, we used a single hidden layer as embedding.
3 DATA: THE 1000 GENOMES PROJECT
The 1000 Genomes project is the first project to sequence the genomes of a large number of people in populations worldwide, yielding the largest public catalog of human genetic variants to date Consortium (2015). This allowed large-scale comparison of DNA sequences from populations, thanks to the presence of genetic variation. Individuals of the 1000 Genomes project are samples taken from 26 populations over the world, which are grouped into 5 geographical regions. Figure 2(a) shows a histogram derived from the 1000 Genomes data, depicting the frequency of individuals per population (ethnicity). Analogously, Figure 2(b) depicts the frequency of individuals per geographical region.
In this dataset, we included 315,345 genetic variants with frequencies of at least 5% in 3,450 individuals sampled worldwide from 26 populations, interrogated using microarray genotyping technology: the Genome-Wide Human SNP Array 6.0 by Affymetrix. The mutated state is established by comparison to the Genome Reference Consortium human genome (build 37). Since individuals have 2 copies of each genomic position, a sampled individual can have 0, 1 or 2 copies of a genetic mutation, hereafter referred to as an individual genotype. We excluded SNPs positioned on the sex
chromosomes and only included SNPs in approximate linkage equilibrium with each other, such that genotypes at neighboring positions are only weakly correlated (r2 < 0.5).
4 EXPERIMENTS
In this section, we describe the model architectures, and report and discuss the obtained results.
4.1 MODEL ARCHITECTURE
We experimented with simple models both in the auxiliary networks and the basic architecture, which yielded very promising results. We designed a basic architecture with 2 hidden layers fol-
lowed by a softmax layer to perform ancestry prediction. We trained this architecture with and without the assistance of the auxiliary network. Similarly, the auxiliary networks were build by stacking a hidden layer on top of one of the feature embeddings described in Section 2.2. In the reported experiments, all hidden layers have 100 units. All models were trained by means of stochastic gradient descent with adaptive learning rate (Tieleman & Hinton, 2012), both for γ = 0 and γ = 10, using dropout, limiting the norm of the weights to 1 and/or applying weight decay to reduce overfitting.2.
4.2 RESULTS
Given the relatively small amount of samples in the 1000 Genomes data, we report results obtained by 5-fold cross validation of the model. We split the data into 5 folds of equal size. A single fold is retained for test, whereas three of the remaining folds are used as training data and the final fold is used as validation data. We repeated the process 5 times (one per fold) and report the means and standard deviations of results on the different test sets.
Table 1 summarizes the results obtained for each model. First, we observe that, for most of Diet Network architectures, training with an reconstruction term in the loss (γ > 0) reduces the misclassification error and provides a lower standard deviation over the folds, suggesting more robustness to variations in the learnt feature embedding.
Training the models end-to-end, with no pre-computed feature embedding, yielded higher misclassification error than simply training the basic model, which could be attributed to the fact that adding the prediction networks makes a difficult, high-dimensional optimization problem even harder. As a general trend, adding pre-computed feature embeddings achieved better performance (lower error), while allowing to significantly reduce the number of free parameters in the fat layers of the model. Among the tested feature embeddings, random projections achieved good results, highlighting the potential of the model when reducing the number of free parameters.
Using the SNP2Vec embedding, trained to exploit the similarities and co-occurences between the SNPs, in conjunction with the Diet Networks framework obtains slightly better results than the model using a random projection. The addition of the reconstruction criterion does not appear to reduce the number of errors made by the model but it does appear to reduce the variance of the results, as observed on the other models.
Despite its simplicity, the per class histogram encoding (when used with a reconstruction criterion) yielded the best results. Note that this encoding is the one with the fewest number of free parameters in the fat layers, with a reduction factor of almost 4000 w.r.t. the analogous basic model (with reconstruction). Figure 3(a) shows the mean results obtained with the histogram embedding. As shown in the figure, when considering the ethnicity, the main misclassifications involve ethnicities likely to display very close genetic proximity, such as British from England and Scotland, and Utah residents with Northern and Western ancestry (likely to be immigrants from England), or Indian Telugu and Sri Lankan Tamil for instance. However, the model achieves almost 100% accuracy when considering the 5 geographical regions.
We also compared the performance of our model to the principal component analysis (PCA) approach, commonly used in the genomics domain, to select subgroups of individuals in order to perform more homogeneous analysis. The number of principal components (PCs) is chosen according to their significance, and usually varies from one dataset to another, being 10 the de facto standard for small datasets. However, in the case of the 1000 Genomes dataset, we could go up to 50 PCs. Therefore, we trained a linear classifier on top of PCA features, considering 10 and 50 PCs, 100 PCs to match the number of feature used in the other experiments, as well as 200 PCs. Using 200 PCs yielded better performance, but going beyond that saturated in terms of misclassification error (see SectionC in Appendix for more details). Adding hidden layers to the classifier didn’t help either (see reported results for several MLP configurations before the linear classifier).
5 CONCLUSION
In this paper, we proposed Diet Networks, a novel network parametrization, which considerably reduces the number of free parameters in the fat layers of a model when the input is very high
2The code to reproduce the experiments can be found here: https://github.com/adri-romsor/DietNetworks
dimensional. We showed how using the parameter prediction networks, yielded better generalization in terms of misclassification error. Notably, when using pre-computed feature embeddings that maximally reduced the number of free parameters, we were able to obtain our best results. We validated our approach on the publicly available 1000 genomes dataset, addressing the relevant task of ancestry prediction based on SNP data. This work demonstrated the potential of deep learning models to tackle domain-specific tasks where there is a mismatch between the number of samples and their high dimensionality.
Given the high accuracy achieved in the ancestry prediction task, we believe that deep learning techniques can improve standard practices in the analysis of human polymorphism data. We expect that these techniques will allow us to tackle the more challenging problem of conducting genetic association studies. Hence, we expect to further develop our method to conduct population-aware analyses of SNP data in disease cohorts. The increased power of deep learning methods to identify the genetic basis of common diseases could lead to better patient risk prediction and will improve our overall understanding of disease etiology.
ACKNOWLEDGMENTS
The authors would like to thank the developers of Theano Theano Development Team (2016) and Lasagne Lasagne (2016). We acknowledge the support of the following agencies for research funding and computing support: Imagia, CIFAR, Canada Research Chairs, Compute Canada and Calcul Québec. J.G.H. is an EPAC/Linacre Junior Research Fellow funded by the Human Frontiers Program (LT-001017/2013-L). Special thanks to Valèria Romero-Soriano, Xavier Grau-Bové and Margaux Luck for their patience sharing genomic biology expertise; as well as to Michal Drozdzal, Caglar Gulcehre and Simon Jégou for useful discussions and support.
A THE 1000 GENOMES PROJECT LEGENDS
A.1 POPULATION ETHNICITY LEGEND
ACB: African Caribbeans in Barbados ASW: Americans of African Ancestry in SW USA BEB: Bengali from Bangladesh CDX: Chinese Dai in Xishuangbanna CEU: Utah Residents (CEPH) with Northern and Western Ancestry CHB: Han Chinese in Bejing CHS: Southern Han Chinese CLM: Colombians from Medellin ESN: Esan in Nigeria FIN: Finnish in Finland GBR: British in England and Scotland GIH: Gujarati Indian from Houston GWD: Gambian in Western Divisions in the Gambia IBS: Iberian Population in Spain ITU: Indian Telugu from the UK JPT: Japanese in Tokyo KHV: Kinh in Ho Chi Minh City LWK: Luhya in Webuye MSL: Mende in Sierra Leone MXL: Mexican Ancestry from Los Angeles PEL: Peruvians from Lima PJL: Punjabi from Lahore PUR: Puerto Ricans STU: Sri Lankan Tamil from the UK TSI: Toscani in Italia YRI: Yoruba in Ibadan
A.2 GEOGPRAHICAL REGION LEGEND
AFR: African AMR: Ad Mixed American EAS: East Asian EUR: European SAS: South Asian
B OBTAINING THE 1000 GENOMES DATA
SNP data for the 1000G dataset was downloaded from ftp://ftp.1000genomes.ebi.ac. uk:21/vol1/ftp/release/20130502/supporting/hd_genotype_chip/
− ALL . wgs . n h g r i c o r i e l l a f f y 6 . 2 0 1 4 0 8 2 5 . g e n o t y p e s h a s p e d . v c f . gz
− a f f y s a m p l e s . 2 0 1 4 1 1 1 8 . p a n e l
Representative commands:
With PLINK v1.90b2n 64-bit https://www.cog-genomics.org/plink2
# c o n v e r t v c f t o p l i n k f o r m a t ( bed ) , and o n l y k e e p i n g common markers ( minor a l l e l e f r e q u e n c y > 0 . 0 5 i n combined sample ) > p l i n k −−v c f $ p a t h /ALL . wgs . n h g r i c o r i e l l a f f y 6 . 2 0 1 4 0 8 2 5 . g e n o t y p e s h a s p e d . v c f . gz −−maf 0 . 0 5 −−o u t $ p a t h / a f f y 6 b i a l l e l i c s n p s m a f 0 0 5 a u t −−not−c h r X Y MT −−make−bed
# produce a pruned s u b s e t o f markers t h a t are i n a p p r o x i m a t e l i n k a g e e q u i l i b r i u m w i t h each o t h e r > p l i n k −−b f i l e $ p a t h / a f f y 6 b i a l l e l i c s n p s m a f 0 0 5 a u t −−indep− p a i r w i s e 50 5 0 . 5 −−o u t $ p a t h / a f f y 6 b i a l l e l i c s n p s m a f 0 0 5 a u t
# e x c l u d e markers t o g e t pruned s u b s e t > p l i n k −−b f i l e $ p a t h / a f f y 6 b i a l l e l i c s n p s m a f 0 0 5 a u t −−e x c l u d e
$ p a t h / a f f y 6 b i a l l e l i c s n p s m a f 0 0 5 a u t . p rune . o u t −−r e c o d e A −− o u t $ p a t h / a f f y 6 b i a l l e l i c s n p s m a f 0 0 5 t h i n n e d a u t A
For information on how to download this pre-processed dataset directly, please email Adriana Romero or Pierre Luc Carrier.
C PCA COMPONENTS AND MISCLASSIFICATION ERROR
In this section, we analyze the influence of increasing the number of PCs used to perform classification. Figure 4 depicts the obtained results when considering 100, 200, 400, 800 and 1000 PCs. As shown in the figure, the best validation error comes with PC200. Further increasing the number of PCs improves the training error but does not generalize well on the validation set.
The paper presents an application of deep learning to genomic SNP data
with a comparison of possible approaches for dealing with the very
high data dimensionality. The approach looks very interesting but the
experiments are too limited to draw firm conclusions about the
strengths of different approaches. The presentation would benefit from
more precise math.
Quality:
The basic idea of the paper is interesting and the applied deep
learning methodology appears reasonable. The experimental evaluation
is rather weak as it only covers a single data set and a very limited
number of cross validation folds. Given the significant variation in
the performances of all the methods, it seems the differences between
the better-performing methods are probably not statistically
significant. More comprehensive empirical validation could clearly
strengthen the paper.
Clarity:
The writing is generally good both in terms of the biology and ML, but
more mathematical rigour would make it easier to understand precisely
what was done. The different architectures are explained on an
intuitive level and might benefit from a clear mathematical
definition. I was ultimately left unsure of what the "raw end2end"
model is - given so few parameters it cannot work on raw 300k
dimensional input but I could not figure out what kind of embedding
was used.
The results in Fig. 3 might be clearer if scaled so that maximum for
each class is 1 to avoid confounding from different numbers of
subjects in different classes. In the text, please use the standard
italics math font for all symbols such as N, N_d, ...
Originality:
The application and the approach appear quite novel.
Significance:
There is clearly strong interest for deep learning in the genomics
area and the paper seeks to address some of the major bottlenecks
here. It is too early to tell whether the specific techniques proposed
in the paper will be the ultimate solution, but at the very least the
paper provides interesting new ideas for others to work on.
Other comments:
I think releasing the code as promised would be a must.