Deep learning models often achieve high accuracy but lack interpretability, making them unsuitable for critical applications such as medical diagnosis, biomolecule design, criminal justice, etc. The Sparse High-order Interaction Model (SHIM) addresses this limitation by providing both transparency and predictive reliability. However, real-world data often contain outliers, which can distort model performance. To overcome this, we propose Huberized-SHIM, an extension of SHIM that integrates Huber loss-based robust regression to mitigate the impact of outliers. We introduce a homotopy-based exact regularization path algorithm and a novel tree-pruning criterion to efficiently manage interaction complexity. Additionally, we incorporate the conformal prediction framework to enhance statistical reliability. Empirical evaluations on synthetic and real-world datasets demonstrate the superior robustness and accuracy of Huberized-SHIM in high-stakes decision-making contexts.

MOTIVATION: Differentiable programming frameworks like PyTorch and JAX revolutionized biological modeling. A foremost merit is that multiple components programmed separately can be put together so that the parameters are jointly optimized. Despite its proven value in agricultural applications, existing breeding simulators are non-differentiable, hindering integration into general deep learning systems. RESULTS: In this paper, we present DiffBreed, a differentiable breeding simulator. Its performance was evaluated in gradient-based optimization of a progeny allocation strategy that maximizes the genetic gain. By utilizing gradient-based optimization, DiffBreed refined progeny allocation strategies, achieving superior genetic gains compared to a non-optimized equal allocation approach. These findings highlight DiffBreed's capacity to properly calculate gradient information through automatic differentiation. With its innovative design, DiffBreed is expected to transform future breeding optimization by integrating into modern deep learning workflows. AVAILABILITY AND IMPLEMENTATION: The proposed automatic differentiation framework was implemented as the Python module "DiffBreed." This module and all the scripts used in this study, including the gradient-based optimization, are available from the "KosukeHamazaki/GORA-PT" repository on GitHub, https://github.com/KosukeHamazaki/GORA-PT. While the simulated datasets in the present study are available from the same repository, the optimized results by PyTorch were not shared in this repository due to the file sizes. Instead, all datasets, including the optimized ones, will be shared in the repository on Zenodo, https://doi.org/10.5281/zenodo.14046522.

Automated high-stake decision-making, such as medical diagnosis, requires models with high interpretability and reliability. We consider the sparse high-order interaction model as an interpretable and reliable model with a good prediction ability. However, finding statistically significant high-order interactions is challenging because of the intrinsically high dimensionality of the combinatorial effects. Another problem in data-driven modeling is the effect of ``cherry-picking" (i.e., selection bias). Our main contribution is extending the recently developed parametric programming approach for selective inference to high-order interaction models. An exhaustive search over the cherry tree (all possible interactions) can be daunting and impractical, even for small-sized problems. We introduced an efficient pruning strategy and demonstrated the computational efficiency and statistical power of the proposed method using both synthetic and real data.

We present a novel nonnegative tensor decomposition method, called Legendre decomposition, which factorizes an input tensor into a multiplicative combination of parameters. Thanks to the well-developed theory of information geometry, the reconstructed tensor is unique and always minimizes the KL divergence from an input tensor. We empirically show that Legendre decomposition can more accurately reconstruct tensors than other nonnegative tensor decomposition methods.

Motivation: Fast and accurate prediction of protein-ligand binding structures is indispensable for structure-based drug design and accurate estimation of binding free energy of drug candidate molecules in drug discovery. Recently, accurate pose prediction methods based on short Molecular Dynamics (MD) simulations, such as MM-PBSA and MM-GBSA, among generated docking poses have been used. Since molecular structures obtained from MD simulation depend on the initial condition, taking the average over different initial conditions leads to better accuracy. Prediction accuracy of protein-ligand binding poses can be improved with multiple runs at different initial velocity. Results: This paper shows that a machine learning method, called Best Arm Identification, can optimally control the number of MD runs for each binding pose. It allows us to identify a correct binding pose with a minimum number of total runs. Our experiment using three proteins and eight inhibitors showed that the computational cost can be reduced substantially without sacrificing accuracy. This method can be applied for controlling all kinds of molecular simulations to obtain best results under restricted computational resources. Availability and implementation: Code and data are available on GitHub at https://github.com/tsudalab/bpbi. Contact: [email protected] or [email protected]. Supplementary information: Supplementary data are available at Bioinformatics online.

Summary: Exhaustive detection of multi-loci markers from genome-wide association study datasets is a computationally challenging problem. This paper presents a massively parallel algorithm for finding all significant combinations of alleles and introduces a software tool termed MP-LAMP that can be easily deployed in a cloud platform, such as Amazon Web Service, as well as in an in-house computer cluster. Multi-loci marker detection is an unbalanced tree search problem that cannot be parallelized by simple tree-splitting using generic parallel programming frameworks, such as Map-Reduce. We employ work stealing and periodic reduce-broadcast to decrease the running time almost linearly to the number of cores. Availability and implementation: MP-LAMP is available at https://github.com/tsudalab/mp-lamp. Supplementary information: Supplementary data are available at Bioinformatics online.

We solve tensor balancing, rescaling an Nth order nonnegative tensor by multiplying N tensors of order N - 1 so that every fiber sums to one. This generalizes a fundamental process of matrix balancing used to compare matrices in a wide range of applications from biology to economics. We present an efficient balancing algorithm with quadratic convergence using Newton’s method and show in numerical experiments that the proposed algorithm is several orders of magnitude faster than existing ones. To theoretically prove the correctness of the algorithm, we model tensors as probability distributions in a statistical manifold and realize tensor balancing as projection onto a submanifold. The key to our algorithm is that the gradient of the manifold, used as a Jacobian matrix in Newton’s method, can be analytically obtained using the Möbius inversion formula, the essential of combinatorial mathematics. Our model is not limited to tensor balancing, but has a wide applicability as it includes various statistical and machine learning models such as weighted DAGs and Boltzmann machines.

Finding statistically significant high-order interactions in predictive modeling is important but challenging task because the possible number of high-order interactions is extremely large (e.g., $> 10^{17}$). In this paper we study feature selection and statistical inference for sparse high-order interaction models. Our main contribution is to extend recently developed selective inference framework for linear models to high-order interaction models by developing a novel algorithm for efficiently characterizing the selection event for the selective inference of high-order interactions. We demonstrate the effectiveness of the proposed algorithm by applying it to an HIV drug response prediction problem.

BACKGROUND: Artificially synthesized RNA molecules provide important ways for creating a variety of novel functional molecules. State-of-the-art RNA inverse folding algorithms can design simple and short RNA sequences of specific GC content, that fold into the target RNA structure. However, their performance is not satisfactory in complicated cases. RESULT: We present a new inverse folding algorithm called MCTS-RNA, which uses Monte Carlo tree search (MCTS), a technique that has shown exceptional performance in Computer Go recently, to represent and discover the essential part of the sequence space. To obtain high accuracy, initial sequences generated by MCTS are further improved by a series of local updates. Our algorithm has an ability to control the GC content precisely and can deal with pseudoknot structures. Using common benchmark datasets for evaluation, MCTS-RNA showed a lot of promise as a standard method of RNA inverse folding. CONCLUSION: MCTS-RNA is available at https://github.com/tsudalab/MCTS-RNA .

We build information geometry for a partially ordered set of variables and define the orthogonal decomposition of information theoretic quantities. The natural connection between information geometry and order theory leads to efficient decomposition algorithms. This generalization of Amari's seminal work on hierarchical decomposition of probability distributions on event combinations enables us to analyze high-order statistical interactions arising in neuroscience, biology, and machine learning.

In this paper we study predictive pattern mining problems where the goal is to construct a predictive model based on a subset of predictive patterns in the database. Our main contribution is to introduce a novel method called safe pattern pruning (SPP) for a class of predictive pattern mining problems. The SPP method allows us to efficiently find a superset of all the predictive patterns in the database that are needed for the optimal predictive model. The advantage of the SPP method over existing boosting-type method is that the former can find the superset by a single search over the database, while the latter requires multiple searches. The SPP method is inspired by recent development of safe feature screening. In order to extend the idea of safe feature screening into predictive pattern mining, we derive a novel pruning rule called safe pattern pruning (SPP) rule that can be used for searching over the tree defined among patterns in the database. The SPP rule has a property that, if a node corresponding to a pattern in the database is pruned out by the SPP rule, then it is guaranteed that all the patterns corresponding to its descendant nodes are never needed for the optimal predictive model. We apply the SPP method to graph mining and item-set mining problems, and demonstrate its computational advantage.

One of the major issues in genome-wide association studies is to solve the missing heritability problem. While considering epistatic interactions among multiple SNPs may contribute to solving this problem, existing software cannot detect statistically significant high-order interactions. We propose software named LAMPLINK, which employs a cutting-edge method to enumerate statistically significant SNP combinations from genome-wide case-control data. LAMPLINK is implemented as a set of additional functions to PLINK, and hence existing procedures with PLINK can be applicable. Applied to the 1000 Genomes Project data, LAMPLINK detected a combination of five SNPs that are statistically significantly accumulated in the Japanese population. AVAILABILITY AND IMPLEMENTATION: LAMPLINK is available at http://a-terada.github.io/lamplink/ CONTACT: [email protected] or [email protected] information: Supplementary data are available at Bioinformatics online.

BACKGROUND: Single-cell RNA sequencing is fast becoming one the standard method for gene expression measurement, providing unique insights into cellular processes. A number of methods, based on general dimensionality reduction techniques, have been suggested to help infer and visualise the underlying structure of cell populations from single-cell expression levels, yet their models generally lack proper biological grounding and struggle at identifying complex differentiation paths. RESULTS: Here we introduce cellTree: an R/Bioconductor package that uses a novel statistical approach, based on document analysis techniques, to produce tree structures outlining the hierarchical relationship between single-cell samples, while identifying latent groups of genes that can provide biological insights. CONCLUSIONS: With cellTree, we provide experimentalists with an easy-to-use tool, based on statistically and biologically-sound algorithms, to efficiently explore and visualise single-cell RNA data. The cellTree package is publicly available in the online Bionconductor repository at: http://bioconductor.org/packages/cellTree/ .

Non-communicable diseases (NCDs) are no longer just a problem for high-income countries, but they are also a problem that affects developing countries. Preventive medicine is definitely the key to combat NCDs; however, the cost of preventive programs is a critical issue affecting the popularization of these medicine programs in developing countries. In this study, we investigate predictive modeling for providing a low-cost preventive medicine program. In our two-year-long field study in Bangladesh, we collected the health checkup results of 15,075 subjects, the data of 6,607 prescriptions, and the follow-up examination results of 2,109 subjects. We address three prediction problems, namely subject risk prediction, drug recommendation, and future risk prediction, by using machine learning techniques; our multiple-classifier approach successfully reduced the costs of health checkups, a multi-task learning method provided accurate recommendation for specific types of drugs, and an active learning method achieved an efficient assignment of healthcare workers for the follow-up care of subjects.

BACKGROUND: Searching for similar compounds in a database is the most important process for in-silico drug screening. Since a query compound is an important starting point for the new drug, a query holder, who is afraid of the query being monitored by the database server, usually downloads all the records in the database and uses them in a closed network. However, a serious dilemma arises when the database holder also wants to output no information except for the search results, and such a dilemma prevents the use of many important data resources. RESULTS: In order to overcome this dilemma, we developed a novel cryptographic protocol that enables database searching while keeping both the query holder's privacy and database holder's privacy. Generally, the application of cryptographic techniques to practical problems is difficult because versatile techniques are computationally expensive while computationally inexpensive techniques can perform only trivial computation tasks. In this study, our protocol is successfully built only from an additive-homomorphic cryptosystem, which allows only addition performed on encrypted values but is computationally efficient compared with versatile techniques such as general purpose multi-party computation. In an experiment searching ChEMBL, which consists of more than 1,200,000 compounds, the proposed method was 36,900 times faster in CPU time and 12,000 times as efficient in communication size compared with general purpose multi-party computation. CONCLUSION: We proposed a novel privacy-preserving protocol for searching chemical compound databases. The proposed method, easily scaling for large-scale databases, may help to accelerate drug discovery research by making full use of unused but valuable data that includes sensitive information.

Three or more transcription factors (TFs) often work together, and the combinatorial regulations are essential in cellular machinery. However, it is impossible to discover statistically significant sets of TF binding motifs due to the necessity of the multiple testing procedure. To improve the sensitivity of widely used Bonferroni correction or its modified methods, such as Holm procedure, Westfall-Young permutation procedure (WY-procedure) has often been applied. However, few studies have used WY-procedure for the discoveries of the combinatorial effects of the motifs because of the extremely large computational time. In this paper, we propose an efficient branch-and-bound algorithm to perform WY-procedure to enumerate statistically significant motif combinations. When we use WY-procedure for the combinatorial regulation discovery, finding the minimum P-value from each permuted dataset consumes an enormous amount of time. We show that a combination that has the possibility to achieve the minimum P-value appears with high frequency over the threshold in dataset. This property enables a frequent itemset mining algorithm to efficiently select the candidates to achieve the minimum P-value. Our demonstrations using yeast and human transcriptome datasets show that the proposed algorithm is orders-of-magnitude faster than WY-procedure, and can practically list statistically significant motif combinations even when any combinations are considered.

Abstract Motivation: Although several methods exist to relate high-dimensional gene expression data to various clinical phenotypes, finding combinations of features in such input remains a challenge, particularly when fitting complex statistical models such as those used for survival studies. Results: Our proposed method builds on existing ‘regularization path-following’ techniques to produce regression models that can extract arbitrarily complex patterns of input features (such as gene combinations) from large-scale data that relate to a known clinical outcome. Through the use of the data’s structure and itemset mining techniques, we are able to avoid combinatorial complexity issues typically encountered with such methods, and our algorithm performs in similar orders of duration as single-variable versions. Applied to data from various clinical studies of cancer patient survival time, our method was able to produce a number of promising gene-interaction candidates whose tumour-related roles appear confirmed by literature. Availability: An R implementation of the algorithm described in this article can be found at https://github.com/david-duverle/regularisation-path-following Contact: [email protected] Supplementary information: Supplementary data are available at Bioinformatics online.

Query autocompletion is an important feature saving users many keystrokes from typing the entire query. In this paper we study the problem of query autocompletion that tolerates errors in users' input using edit distance constraints. Previous approaches index data strings in a trie, and continuously maintain all the prefixes of data strings whose edit distance from the query are within the threshold. The major inherent problem is that the number of such prefixes is huge for the first few characters of the query and is exponential in the alphabet size. This results in slow query response even if the entire query approximately matches only few prefixes. In this paper, we propose a novel neighborhood generation-based algorithm, IncNGTrie, which can achieve up to two orders of magnitude speedup over existing methods for the error-tolerant query autocompletion problem. Our proposed algorithm only maintains a small set of active nodes, thus saving both space and time to process the query. We also study efficient duplicate removal which is a core problem in fetching query answers. In addition, we propose optimization techniques to reduce our index size, as well as discussions on several extensions to our method. The efficiency of our method is demonstrated against existing methods through extensive experiments on real datasets.

Similarity search in databases of labeled graphs is a fundamental task in managing graph data such as XML, chemical compounds and social networks. Typically, a graph is decomposed to a set of substructures (e.g., paths, trees and subgraphs) and a similarity measure is defined via the number of common substructures. Using the representation, graphs can be stored in a document database by regarding graphs as documents and substructures as words. A graph similarity query then translates to a semi-conjunctive query that retrieves graphs sharing at least k substructures in common with the query graph. We argue that this kind of query cannot be solved efficiently by conventional inverted indexes, and develop a novel recursive search algorithm on wavelet trees (Grossi et al., SODA'03). Unlike gIndex, it does not require frequent subgraph mining for indexing. In experiments, our method was successfully applied to 25 million chemical compounds.

MOTIVATION: Recent progress in DNA sequencing technologies calls for fast and accurate algorithms that can evaluate sequence similarity for a huge amount of short reads. Searching similar pairs from a string pool is a fundamental process of de novo genome assembly, genome-wide alignment and other important analyses. RESULTS: In this study, we designed and implemented an exact algorithm SlideSort that finds all similar pairs from a string pool in terms of edit distance. Using an efficient pattern growth algorithm, SlideSort discovers chains of common k-mers to narrow down the search. Compared to existing methods based on single k-mers, our method is more effective in reducing the number of edit distance calculations. In comparison to backtracking methods such as BWA, our method is much faster in finding remote matches, scaling easily to tens of millions of sequences. Our software has an additional function of single link clustering, which is useful in summarizing short reads for further processing. AVAILABILITY: Executable binary files and C++ libraries are available at http://www.cbrc.jp/~shimizu/slidesort/ for Linux and Windows.

The analysis of n -ary relations receives attention in many different fields, for instance biology, web mining, and social studies. In the basic setting, there are n sets of instances, and each observation associates n instances, one from each set. A common approach to explore these n -way data is the search for n -set patterns, the n -way equivalent of itemsets. More precisely, an n -set pattern consists of specific subsets of the n instance sets such that all possible associations between the corresponding instances are observed in the data. In contrast, traditional itemset mining approaches consider only two-way data, namely items versus transactions. The n -set patterns provide a higher-level view of the data, revealing associative relationships between groups of instances. Here, we generalize this approach in two respects. First, we tolerate missing observations to a certain degree, that means we are also interested in n -sets where most (although not all) of the possible associations have been recorded in the data. Second, we take association weights into account. In fact, we propose a method to enumerate all n -sets that satisfy a minimum threshold with respect to the average association weight. Technically, we solve the enumeration task using a reverse search strategy, which allows for effective pruning of the search space. In addition, our algorithm provides a ranking of the solutions and can consider further constraints. We show experimental results on artificial and real-world datasets from different domains.

BACKGROUND: Understanding of secondary metabolic pathway in plant is essential for finding druggable candidate enzymes. However, there are many enzymes whose functions are not yet discovered in organism-specific metabolic pathways. Towards identifying the functions of those enzymes, assignment of EC numbers to the enzymatic reactions they catalyze plays a key role, since EC numbers represent the categorization of enzymes on one hand, and the categorization of enzymatic reactions on the other hand. RESULTS: We propose reaction graph kernels for automatically assigning EC numbers to unknown enzymatic reactions in a metabolic network. Reaction graph kernels compute similarity between two chemical reactions considering the similarity of chemical compounds in reaction and their relationships. In computational experiments based on the KEGG/REACTION database, our method successfully predicted the first three digits of the EC number with 83% accuracy. We also exhaustively predicted missing EC numbers in plant's secondary metabolism pathway. The prediction results of reaction graph kernels on 36 unknown enzymatic reactions are compared with an expert's knowledge. Using the same data for evaluation, we compared our method with E-zyme, and showed its ability to assign more number of accurate EC numbers. CONCLUSION: Reaction graph kernels are a new metric for comparing enzymatic reactions.

Several key problems in machine learning, such as feature selection and active learning, can be formulated as submodular set function maximization. We present herein a novel algorithm for maximizing a submodular set function under a cardinality constraint --- the algorithm is based on a cutting-plane method and is implemented as an iterative small-scale binary-integer linear programming procedure. It is well known that this problem is NP-hard, and the approximation factor achieved by the greedy algorithm is the theoretical limit for polynomial time. As for (non-polynomial time) exact algorithms that perform reasonably in practice, there has been very little in the literature although the problem is quite important for many applications. Our algorithm is guaranteed to find the exact solution in finite iterations, and it converges fast in practice due to the efficiency of the cutting-plane mechanism. Moreover, we also provide a method that produces successively decreasing upper-bounds of the optimal solution, while our algorithm provides successively increasing lower-bounds. Thus, the accuracy of the current solution can be estimated at any point, and the algorithm can be stopped early once a desired degree of tolerance is met. We evaluate our algorithm on sensor placement and feature selection applications showing good performance.

Many real-world applications with graph data require the solution of a given regression task as well as the identification of the subgraphs which are relevant for the task. In these cases graphs are commonly represented as high dimensional binary vectors of indicators of subgraphs. However, since the dimensionality of such indicator vectors can be high even for small datasets, traditional regression algorithms become intractable and past approaches used to preselect a feasible subset of subgraphs. A different approach was recently proposed by a Lasso-type method where the objective function optimization with a large number of variables is reformulated as a dual mathematical programming problem with a small number of variables but a large number of constraints. The dual problem is then solved by column generation, where the subgraphs corresponding to the most violated constraints are found by weighted subgraph mining. This paper proposes an extension of this method to a Bayesian approach in which the regression parameters are considered as random variables and integrated out from the model likelihood, thus providing a posterior distribution on the target variable as opposed to a point estimate. We focus on a linear regression model with a Gaussian prior distribution on the parameters. We evaluate our approach on several molecular graph datasets and analyze whether the uncertainty in the target estimate given by the target posterior distribution variance can be used to improve model performance and therefore provides useful additional information.

We propose Link Propagation as a new semi-supervised learning method for link prediction problems, where the task is to predict unknown parts of the network structure by using auxiliary information such as node similarities. Since the proposed method can fill in missing parts of tensors, it is applicable to multi-relational domains, allowing us to handle multiple types of links simultaneously. We also give a novel efficient algorithm for Link Propagation based on an accelerated conjugate gradient method.

SUMMARY: The DICS database is a dynamic web repository of computationally predicted functional modules from the human protein-protein interaction network. It provides references to the CORUM, DrugBank, KEGG and Reactome pathway databases. DICS can be accessed for retrieving sets of overlapping modules and protein complexes that are significantly enriched in a gene list, thereby providing valuable information about the functional context. AVAILABILITY: Supplementary information on datasets and methods is available on the web server http://mips.gsf.de/proj/dics.

MOTIVATION: Modern systems biology aims at understanding how the different molecular components of a biological cell interact. Often, cellular functions are performed by complexes consisting of many different proteins. The composition of these complexes may change according to the cellular environment, and one protein may be involved in several different processes. The automatic discovery of functional complexes from protein interaction data is challenging. While previous approaches use approximations to extract dense modules, our approach exactly solves the problem of dense module enumeration. Furthermore, constraints from additional information sources such as gene expression and phenotype data can be integrated, so we can systematically mine for dense modules with interesting profiles. RESULTS: Given a weighted protein interaction network, our method discovers all protein sets that satisfy a user-defined minimum density threshold. We employ a reverse search strategy, which allows us to exploit the density criterion in an efficient way. Our experiments show that the novel approach is feasible and produces biologically meaningful results. In comparative validation studies using yeast data, the method achieved the best overall prediction performance with respect to confirmed complexes. Moreover, by enhancing the yeast network with phenotypic and phylogenetic profiles and the human network with tissue-specific expression data, we identified condition-dependent complex variants. AVAILABILITY: A C++ implementation of the algorithm is available at http://www.kyb.tuebingen.mpg.de/~georgii/dme.html. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

MOTIVATION: The existing supervised methods for biological network inference work on each of the networks individually based only on intra-species information such as gene expression data. We believe that it will be more effective to use genomic data and cross-species evolutionary information from different species simultaneously, rather than to use the genomic data alone. RESULTS: We created a new semi-supervised learning method called Link Propagation for inferring biological networks of multiple species based on genome-wide data and evolutionary information. The new method was applied to simultaneous reconstruction of three metabolic networks of Caenorhabditis elegans, Helicobacter pylori and Saccharomyces cerevisiae, based on gene expression similarities and amino acid sequence similarities. The experimental results proved that the new simultaneous network inference method consistently improves the predictive performance over the individual network inferences, and it also outperforms in accuracy and speed other established methods such as the pairwise support vector machine. AVAILABILITY: The software and data are available at http://cbio.ensmp.fr/~yyamanishi/LinkPropagation/.

MOTIVATION: We address the issue of finding a three-way gene interaction, i.e. two interacting genes in expression under the genotypes of another gene, given a dataset in which expressions and genotypes are measured at once for each individual. This issue can be a general, switching mechanism in expression of two genes, being controlled by categories of another gene, and finding this type of interaction can be a key to elucidating complex biological systems. The most suitable method for this issue is likelihood ratio test using logistic regressions, which we call interaction test, but a serious problem of this test is computational intractability at a genome-wide level. RESULTS: We developed a fast method for this issue which improves the speed of interaction test by around 10 times for any size of datasets, keeping highly interacting genes with an accuracy of approximately 85%. We applied our method to approximately 3 x 10(8) three-way combinations generated from a dataset on human brain samples and detected three-way gene interactions with small P-values. To check the reliability of our results, we first conducted permutations by which we can show that the obtained P-values are significantly smaller than those obtained from permuted null examples. We then used GEO (Gene Expression Omnibus) to generate gene expression datasets with binary classes to confirm the detected three-way interactions by using these datasets and interaction tests. The result showed us some datasets with significantly small P-values, strongly supporting the reliability of the detected three-way interactions. AVAILABILITY: Software is available from http://www.bic.kyoto-u.ac.jp/pathway/kayano/bioinfo_three-way.html CONTACT: [email protected] SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

Graph mining methods enumerate frequently appearing subgraph patterns, which can be used as features for subsequent classification or regression. However, frequent patterns are not necessarily informative for the given learning problem. We propose a mathematical programming boosting method (gBoost) that progressively collects informative patterns. Compared to AdaBoost, gBoost can build the prediction rule with fewer iterations. To apply the boosting method to graph data, a branch-and-bound pattern search algorithm is developed based on the DFS code tree. The constructed search space is reused in later iterations to minimize the computation time. Our method can learn more efficiently than the simpler method based on frequent substructure mining, because the output labels are used as an extra information source for pruning the search space. Furthermore, by engineering the mathematical program, a wide range of machine learning problems can be solved without modifying the pattern search algorithm.

Discovery of knowledge from geometric graph databases is of particular importance in chemistry and biology, because chemical compounds and proteins are represented as graphs with 3D geometric coordinates. In such applications, scientists are not interested in the statistics of the whole database. Instead they need information about a novel drug candidate or protein at hand, represented as a query graph. We propose a polynomial-delay algorithm for geometric frequent subgraph retrieval. It enumerates all subgraphs of a single given query graph which are frequent geometric epsi-subgraphs under the entire class of rigid geometric transformations in a database. By using geometric epsi-subgraphs, we achieve tolerance against variations in geometry. We compare the proposed algorithm to gSpan on chemical compound data, and we show that for a given minimum support the total number of frequent patterns is substantially limited by requiring geometric matching. Although the computation time per pattern is larger than for non-geometric graph mining, the total time is within a reasonable level even for small minimum support.

Graph mining methods enumerate frequent subgraphs efficiently, but they are not necessarily good features for machine learning due to high correlation among features. Thus it makes sense to perform principal component analysis to reduce the dimensionality and create decorrelated features. We present a novel iterative mining algorithm that captures informative patterns corresponding to major entries of top principal components. It repeatedly calls weighted substructure mining where example weights are updated in each iteration. The Lanczos algorithm, a standard algorithm of eigen decomposition, is employed to update the weights. In experiments, our patterns are shown to approximate the principal components obtained by frequent mining.

Attributed graphs are increasingly more common in many application domains such as chemistry, biology and text processing. A central issue in graph mining is how to collect informative subgraph patterns for a given learning task. We propose an iterative mining method based on partial least squares regression (PLS). To apply PLS to graph data, a sparse version of PLS is developed first and then it is combined with a weighted pattern mining algorithm. The mining algorithm is iteratively called with different weight vectors, creating one latent component per one mining call. Our method, graph PLS, is efficient and easy to implement, because the weight vector is updated with elementary matrix calculations. In experiments, our graph PLS algorithm showed competitive prediction accuracies in many chemical datasets and its efficiency was significantly superior to graph boosting (gBoost) and the naive method based on frequent graph mining.

Graph data such as chemical compounds and XML documents are getting more common in many application domains. A main difficulty of graph data processing lies in the intrinsic high dimensionality of graphs, namely, when a graph is represented as a binary feature vector of indicators of all possible subgraph patterns, the dimensionality gets too large for usual statistical methods. We propose a nonparametric Bayesian method for clustering graphs and selecting salient patterns at the same time. Variational inference is adopted here, because sampling is not applicable due to extremely high dimensionality. The feature set minimizing the free energy is efficiently collected with the DFS code tree, where the generation of useless subgraphs is suppressed by a tree pruning condition. In experiments, our method is compared with a simpler approach based on frequent subgraph mining, and graph kernels.

In Web-related applications of image categorization, it is desirable to derive an interpretable classification rule with high accuracy. Using the bag-of-words representation and the linear support vector machine, one can partly fulfill the goal, but the accuracy of linear classifiers is not high and the obtained features are not informative for users. We propose to combine item set mining and large margin classifiers to select features from the power set of all visual words. Our resulting classification rule is easier to browse and simpler to understand, because each feature has richer information. As a next step, each image is represented as a graph where nodes correspond to local image features and edges encode geometric relations between features. Combining graph mining and boosting, we can obtain a classification rule based on subgraph features that contain more information than the set features. We evaluate our algorithm in a web-retrieval ranking task where the goal is to reject outliers from a set of images returned for a keyword query. Furthermore, it is evaluated on the supervised classification tasks with the challenging VOC2005 data set. Our approach yields excellent accuracy in the unsupervised ranking task compared to a recently proposed probabilistic model and competitive results in the supervised classification task.

Recent approaches to action classification in videos have used sparse spatio-temporal words encoding local appearance around interesting movements. Most of these approaches use a histogram representation, discarding the temporal order among features. But this ordering information can contain important information about the action itself e.g. consider the sport disciplines of hurdle race and long jump, where the global temporal order of motions (running, jumping) is important to discriminate between the two. In this work we propose to use a sequential representation which retains this temporal order. Further, we introduce Discriminative Subsequence Mining to find optimal discriminative subsequence patterns. In combination with the LPBoost classifier, this amounts to simultaneously learning a classification function and performing feature selection in the space of all possible feature sequences. The resulting classifier linearly combines a small number of interpretable decision functions, each checking for the presence of a single discriminative pattern. The classifier is benchmarked on the KTH action classification data set and outperforms the best known results in the literature.

Graph data such as chemical compounds and XML documents are getting more common in many application domains. A main difficulty of graph data processing lies in the intrinsic high dimensionality of graphs, namely, when a graph is represented as a binary feature vector of indicators of all possible subgraph patterns, the dimensionality gets too large for usual statistical methods. We propose an efficient method to select a small number of salient patterns by regularization path tracking. The generation of useless patterns is minimized by progressive extension of the search space. In experiments, it is shown that our technique is considerably more efficient than a simpler approach based on frequent substructure mining.

We propose an efficient algorithm for principal component analysis (PCA) that is applicable when only the inner product with a given vector is needed. We show that Krylov subspace learning works well both in matrix compression and implicit calculation of the inner product by taking full advantage of the arbitrariness of the seed vector. We apply our algorithm to a PCA-based change-point detection algorithm, and show that it results in about 50 times improvement in computational time.

MOTIVATION: Human immunodeficiency virus type 1 (HIV-1) evolves in human body, and its exposure to a drug often causes mutations that enhance the resistance against the drug. To design an effective pharmacotherapy for an individual patient, it is important to accurately predict the drug resistance based on genotype data. Notably, the resistance is not just the simple sum of the effects of all mutations. Structural biological studies suggest that the association of mutations is crucial: even if mutations A or B alone do not affect the resistance, a significant change might happen when the two mutations occur together. Linear regression methods cannot take the associations into account, while decision tree methods can reveal only limited associations. Kernel methods and neural networks implicitly use all possible associations for prediction, but cannot select salient associations explicitly. RESULTS: Our method, itemset boosting, performs linear regression in the complete space of power sets of mutations. It implements a forward feature selection procedure where, in each iteration, one mutation combination is found by an efficient branch-and-bound search. This method uses all possible combinations, and salient associations are explicitly shown. In experiments, our method worked particularly well for predicting the resistance of nucleotide reverse transcriptase inhibitors (NRTIs). Furthermore, it successfully recovered many mutation associations known in biological literature. AVAILABILITY: http://www.kyb.mpg.de/bs/people/hiroto/iboost/. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

The field of computational biology has seen dramatic growth over the past few years, both in terms of available data, scientific questions and challenges for learning and inference. These new types of scientific and clinical problems require the development of novel supervised and unsupervised learning approaches. In particular, the field is characterized by a diversity of heterogeneous data. The human genome sequence is accompanied by real-valued gene expression data, functional annotation of genes, genotyping information, a graph of interacting proteins, a set of equations describing the dynamics of a system, localization of proteins in a cell, a phylogenetic tree relating species, natural language text in the form of papers describing experiments, partial models that provide priors, and numerous other data sources. This supplementary issue consists of seven peer-reviewed papers based on the NIPS Workshop on New Problems and Methods in Computational Biology held at Whistler, British Columbia, Canada on December 8, 2006. The Neural Information Processing Systems Conference is the premier scientific meeting on neural computation, with session topics spanning artificial intelligence, learning theory, neuroscience, etc. The goal of this workshop was to present emerging problems and machine learning techniques in computational biology, with a particular emphasis on methods for computational learning from heterogeneous data. We received 37 extended abstract submissions, from which 13 were selected for oral presentation. The current supplement contains seven papers based on a subset of the 13 extended abstracts. Submitted manuscripts were rigorously reviewed by at least two referees. The quality of each paper was evaluated with respect to its contribution to biology as well as the novelty of the machine learning methods employed.

We analysed human interaction data from MINT, Intact, HPRD, and DIP in the context of tissue-specific gene expression data in human provided by Su et al . [ 3 ]. We discretized the expression information into binary states (expressed versus not expressed) and searched for densely connected modules where all proteins are expressed in at least 3 tissues and all proteins are not expressed in at least 10 tissues. To deal with the fact that protein interaction data contain a high number of false positives, we computed reliability scores for each experimental source. Similarly to the work by Jansen et al . [ 4 ], we used for that purpose a gold standard set of known interactions as well as a gold standard set of false interactions and calculated the likelihood ratio, which was used to assign edge weights to the interaction graph. The density of a module is defined as the sum of the edge weights inside the module divided by the maximal possible weight sum for a module of that size. Setting the minimum density threshold to 35% and removing modules that are totally contained in other modules, we obtained a set of 949 differentially expressed modules. They were ranked in descending order according to the average weight per node (see [ 5 ]), so larger and denser modules appear first. On the one hand, we discovered known complexes and modules that link strongly cooperating complexes like MCM and ORC. On the other hand, we found extensions of known complexes that confirm hypothetical functional annotation in Uniprot as well as modules which are not contained in the manually curated set of known complexes, but share the same functional annotation. Finally, some modules are candidates for further biological investigation, containing proteins with unknown functional relationships. We developed a general method for exhaustive dense module extraction from networks. Remarkably, it allows to determine exact P-values for the predicted modules without having to rely on any network model and can easily integrate information from different heterogeneous data sources.

Graph data is getting increasingly popular in, e.g., bioinformatics and text processing. A main difficulty of graph data processing lies in the intrinsic high dimensionality of graphs, namely, when a graph is represented as a binary feature vector of indicators of all possible subgraphs, the dimensionality gets too large for usual statistical methods. We propose an efficient method for learning a binomial mixture model in this feature space. Combining the l1 regularizer and the data structure called DFS code tree, the MAP estimate of non-zero parameters are computed efficiently by means of the EM algorithm. Our method is applied to the clustering of RNA graphs, and is compared favorably with graph kernels and the spectral graph distance.

Abstract Motivation: In detection of non-coding RNAs, it is often necessary to identify the secondary structure motifs from a set of putative RNA sequences. Most of the existing algorithms aim to provide the best motif or few good motifs, but biologists often need to inspect all the possible motifs thoroughly. Results: Our method RNAmine employs a graph theoretic representation of RNA sequences and detects all the possible motifs exhaustively using a graph mining algorithm. The motif detection problem boils down to finding frequently appearing patterns in a set of directed and labeled graphs. In the tasks of common secondary structure prediction and local motif detection from long sequences, our method performed favorably both in accuracy and in efficiency with the state-of-the-art methods such as CMFinder. Availability: The software is available upon request. Contact: [email protected] Supplementary information: Visit the following URL for Supplementary information, software availability and the information about the web server:

MOTIVATION: The functions of non-coding RNAs are strongly related to their secondary structures, but it is known that a secondary structure prediction of a single sequence is not reliable. Therefore, we have to collect similar RNA sequences with a common secondary structure for the analyses of a new non-coding RNA without knowing the exact secondary structure itself. Therefore, the sequence comparison in searching similar RNAs should consider not only their sequence similarities but also their potential secondary structures. Sankoff's algorithm predicts the common secondary structures of the sequences, but it is computationally too expensive to apply to large-scale analyses. Because we often want to compare a large number of cDNA sequences or to search similar RNAs in the whole genome sequences, much faster algorithms are required. RESULTS: We propose a new method of comparing RNA sequences based on the structural alignments of the fixed-length fragments of the stem candidates. The implemented software, SCARNA (Stem Candidate Aligner for RNAs), is fast enough to apply to the long sequences in the large-scale analyses. The accuracy of the alignments is better or comparable with the much slower existing algorithms. AVAILABILITY: The web server of SCARNA with graphical structural alignment viewer is available at http://www.scarna.org/.

BACKGROUND: Prediction of human cell response to anti-cancer drugs (compounds) from microarray data is a challenging problem, due to the noise properties of microarrays as well as the high variance of living cell responses to drugs. Hence there is a strong need for more practical and robust methods than standard methods for real-value prediction. RESULTS: We devised an extended version of the off-subspace noise-reduction (de-noising) method to incorporate heterogeneous network data such as sequence similarity or protein-protein interactions into a single framework. Using that method, we first de-noise the gene expression data for training and test data and also the drug-response data for training data. Then we predict the unknown responses of each drug from the de-noised input data. For ascertaining whether de-noising improves prediction or not, we carry out 12-fold cross-validation for assessment of the prediction performance. We use the Pearson's correlation coefficient between the true and predicted response values as the prediction performance. De-noising improves the prediction performance for 65% of drugs. Furthermore, we found that this noise reduction method is robust and effective even when a large amount of artificial noise is added to the input data. CONCLUSION: We found that our extended off-subspace noise-reduction method combining heterogeneous biological data is successful and quite useful to improve prediction of human cell cancer drug responses from microarray data.

In the information regularization framework by Corduneanu and Jaakkola (2005), the distributions of labels are propagated on a hypergraph for semi-supervised learning. The learning is efficiently done by a Blahut-Arimoto-like two step algorithm, but, unfortunately, one of the steps cannot be solved in a closed form. In this paper, we propose a dual version of information regularization, which is considered as more natural in terms of information geometry. Our learning algorithm has two steps, each of which can be solved in a closed form. Also it can be naturally applied to exponential family distributions such as Gaussians. In experiments, our algorithm is applied to protein classification based on a metabolic network and known functional categories.

MOTIVATION: Inferring networks of proteins from biological data is a central issue of computational biology. Most network inference methods, including Bayesian networks, take unsupervised approaches in which the network is totally unknown in the beginning, and all the edges have to be predicted. A more realistic supervised framework, proposed recently, assumes that a substantial part of the network is known. We propose a new kernel-based method for supervised graph inference based on multiple types of biological datasets such as gene expression, phylogenetic profiles and amino acid sequences. Notably, our method assigns a weight to each type of dataset and thereby selects informative ones. Data selection is useful for reducing data collection costs. For example, when a similar network inference problem must be solved for other organisms, the dataset excluded by our algorithm need not be collected. RESULTS: First, we formulate supervised network inference as a kernel matrix completion problem, where the inference of edges boils down to estimation of missing entries of a kernel matrix. Then, an expectation-maximization algorithm is proposed to simultaneously infer the missing entries of the kernel matrix and the weights of multiple datasets. By introducing the weights, we can integrate multiple datasets selectively and thereby exclude irrelevant and noisy datasets. Our approach is favorably tested in two biological networks: a metabolic network and a protein interaction network. AVAILABILITY: Software is available on request.

We address the problem of learning a symmetric positive definite matrix. The central issue is to design parameter updates that preserve positive definiteness. Our updates are motivated with the von Neumann divergence. Rather than treating the most general case, we focus on two key applications that exemplify our methods: on-line learning with a simple square loss, and finding a symmetric positive definite matrix subject to linear constraints. The updates generalize the exponentiated gradient (EG) update and AdaBoost, respectively: the parameter is now a symmetric positive definite matrix of trace one instead of a probability vector (which in this context is a diagonal positive definite matrix with trace one). The generalized updates use matrix logarithms and exponentials to preserve positive definiteness. Most importantly, we show how the derivation and the analyses of the original EG update and AdaBoost generalize to the non-diagonal case. We apply the resulting matrix exponentiated gradient (MEG) update and DefiniteBoost to the problem of learning a kernel matrix from distance measurements.

One way of image denoising is to project a noisy image to the subspace of admissible images derived, for instance, by PCA. However, a major drawback of this method is that all pixels are updated by the projection, even when only a few pixels are corrupted by noise or occlusion. We propose a new method to identify the noisy pixels by l1-norm penalization and to update the identified pixels only. The identification and updating of noisy pixels are formulated as one linear program which can be efficiently solved. In particular, one can apply the upsilon trick to directly specify the fraction of pixels to be reconstructed. Moreover, we extend the linear program to be able to exploit prior knowledge that occlusions often appear in contiguous blocks (e.g., sunglasses on faces). The basic idea is to penalize boundary points and interior points of the occluded area differently. We are also able to show the upsilon property for this extended LP leading to a method which is easy to use. Experimental results demonstrate the power of our approach.

We address the problem of learning a symmetric positive definite matrix. The central issue is to design parameter updates that preserve positive definiteness. Our updates are motivated with the von Neumann diver- gence. Rather than treating the most general case, we focus on two key applications that exemplify our methods: On-line learning with a simple square loss and finding a symmetric positive definite matrix subject to symmetric linear constraints. The updates generalize the Exponentiated Gradient (EG) update and AdaBoost, respectively: the parameter is now a symmetric positive definite matrix of trace one instead of a probability vector (which in this context is a diagonal positive definite matrix with trace one). The generalized updates use matrix logarithms and exponen- tials to preserve positive definiteness. Most importantly, we show how the analysis of each algorithm generalizes to the non-diagonal case. We apply both new algorithms, called the Matrix Exponentiated Gradient (MEG) update and DefiniteBoost, to learn a kernel matrix from distance measurements. 1 Introduction Most learning algorithms have been developed to learn a vector of parameters from data. However, an increasing number of papers are now dealing with more structured parame- ters. More specifically, when learning a similarity or a distance function among objects, the parameters are defined as a symmetric positive definite matrix that serves as a kernel (e.g. [14, 11, 13]). Learning is typically formulated as a parameter updating procedure to optimize a loss function. The gradient descent update [6] is one of the most commonly used algorithms, but it is not appropriate when the parameters form a positive definite matrix, because the updated parameter is not necessarily positive definite. Xing et al. [14] solved this problem by always correcting the updated matrix to be positive. However no bound has been proven for this update-and-correction approach. In this paper, we introduce the Matrix Exponentiated Gradient update which works as follows: First, the matrix logarithm of the current parameter matrix is computed. Then a step is taken in the direction of the steepest descent. Finally, the parameter matrix is updated to the exponential of the modified log-matrix. Our update preserves symmetry and positive definiteness because the matrix exponential maps any symmetric matrix to a positive definite matrix. Bregman divergences play a central role in the motivation and the analysis of on-line learn- ing algorithms [5]. A learning problem is essentially defined by a loss function, and a di- vergence that measures the discrepancy between parameters. More precisely, the updates are motivated by minimizing the sum of the loss function and the Bregman divergence, where the loss function is multiplied by a positive learning rate. Different divergences lead to radically different updates [6]. For example, the gradient descent is derived from the squared Euclidean distance, and the exponentiated gradient from the Kullback-Leibler di- vergence. We use the von Neumann divergence (also called quantum relative entropy) for measuring the discrepancy between two positive definite matrices [8]. We derive a new Matrix Exponentiated Gradient update from this divergence (which is a Bregman diver- gence for positive definite matrices). Finally we prove relative loss bounds using the von Neumann divergence as a measure of progress. Also the following related key problem has received a lot of attention recently [14, 11, 13]: Find a symmetric positive definite matrix that satisfies a number of symmetric linear inequality constraints. The new DefiniteBoost algorithm greedily chooses the most violated constraint and performs an approximated Bregman projection. In the diagonal case, we recover AdaBoost [9]. We also show how the convergence proof of AdaBoost generalizes to the non-diagonal case. 2 von Neumann Divergence or Quantum Relative Entropy If F is a real convex differentiable function on the parameter domain (symmetric d d positive definite matrices) and f (W) := F(W), then the Bregman divergence between two parameters W and W is defined as F(W, W) = F(W) - F(W) - tr[(W - W)f(W)]. When choosing F(W) = tr(W log W -W), then f(W) = log W and the corresponding Bregman divergence becomes the von Neumann divergence [8]: F(W, W) = tr(W log W - W log W - W + W). (1) In this paper, we are primarily interested in the normalized case (when tr(W) = 1). In this case, the positive symmetric definite matrices are related to density matrices commonly used in Statistical Physics and the divergence simplifies to F(W, W) = tr(W log W - W log W). If W = i ivivi is our notation for the eigenvalue decomposition, then we can rewrite the normalized divergence as ~ ~ F(W, W) = i ln ~ i + i ln j(~ vi vj )2. i i,j So this divergence quantifies the difference in the eigenvalues as well as the eigenvectors. 3 On-line Learning In this section, we present a natural extension of the Exponentiated Gradient (EG) up- date [6] to an update for symmetric positive definite matrices. At the t-th trial, the algorithm receives a symmetric instance matrix Xt Rdd. It then produces a prediction ^ yt = tr(WtXt) based on the algorithm's current symmetric positive definite parameter matrix Wt. Finally it incurs for instance1 a quadratic loss (^ yt - yt)2, 1For the sake of simplicity, we use the simple quadratic loss: L W t (W) = (tr(Xt ) - yt)2. For the general update, the gradient Lt(Wt) is exponentiated in the update (4) and this gradient must be symmetric. Following [5], more general loss functions (based on Bregman divergences) are amenable to our techniques. and updates its parameter matrix Wt. In the update we aim to solve the following problem: Wt+1 = argmin W F(W, Wt) + (tr(WXt) - yt)2 , (2) where the convex function F defines the Bregman divergence. Setting the derivative with respect to W to zero, we have f (Wt+1) - f(Wt) + [(tr(Wt+1Xt) - yt)2] = 0. (3) The update rule is derived by solving (3) with respect to Wt+1, but it is not solvable in closed form. A common way to avoid this problem is to approximate tr(Wt+1Xt) by tr(WtXt) [5]. Then, we have the following update: Wt+1 = f -1(f (Wt) - 2(^yt - yt)Xt). In our case, F(W) = tr(W log W - W) and thus f(W) = log W and f-1(W) = exp W. We also augment (2) with the constraint tr(W) = 1, leading to the following Matrix Exponential Gradient (MEG) Update: 1 Wt+1 = exp(log W Z t - 2(^yt - yt)Xt), (4) t where the normalization factor Zt is tr[exp(log Wt - 2(^yt - yt)Xt)]. Note that in the above update, the exponent log Wt - 2(^yt - yt)Xt is an arbitrary symmetric matrix and the matrix exponential converts this matrix back into a symmetric positive definite matrix. A numerically stable version of the MEG update is given in Section 3.2. 3.1 Relative Loss Bounds We now begin with the definitions needed for the relative loss bounds. Let S = (X1, y1), . . . , (XT , yT ) denote a sequence of examples, where the instance matrices Xt Rdd are symmetric and the labels yt R. For any symmetric positive semi-definite ma- trix U with tr(U) = 1, define its total loss as LU(S) = T (tr(UX t=1 t ) - yt)2. The total loss of the on-line algorithm is LMEG(S) = T (tr(W t=1 tXt) - yt)2. We prove a bound on the relative loss LMEG(S) -LU(S) that holds for any U. The proof generalizes a sim- ilar bound for the Exponentiated Gradient update (Lemmas 5.8 and 5.9 of [6]). The relative loss bound is derived in two steps: Lemma 3.1 bounds the relative loss for an individual trial and Lemma 3.2 for a whole sequence (Proofs are given in the full paper). Lemma 3.1 Let Wt be any symmetric positive definite matrix. Let Xt be any symmetric matrix whose smallest and largest eigenvalues satisfy max - min r. Assume Wt+1 is produced from Wt by the MEG update and let U be any symmetric positive semi-definite matrix. Then for any constants a and b such that 0 < a 2b/(2 + r2b) and any learning rate = 2b/(2 + r2b), we have a(yt - tr(WtXt))2 - b(yt - tr(UXt))2 (U,Wt) - (U,Wt+1) (5) In the proof, we use the Golden-Thompson inequality [3], i.e., tr[exp(A + B)] tr[exp(A) exp(B)] for symmetric matrices A and B. We also needed to prove the fol- lowing generalization of Jensen's inequality to matrices: exp(1A + 2(I - A)) exp(1)A + exp(2)(I - A) for finite 1,2 R and any symmetric matrix A with 0 < A I. These two key inequalities will also be essential for the analysis of Definite- Boost in the next section. Lemma 3.2 Let W1 and U be arbitrary symmetric positive definite initial and comparison matrices, respectively. Then for any c such that = 2c/(r2(2 + c)), c 1 1 LMEG(S) 1 + LU(S) + + r2(U, W 2 2 c 1). (6) Proof For the maximum tightness of (5), a should be chosen as a = = 2b/(2 + r2b). Let b = c/r2, and thus a = 2c/(r2(2 + c)). Then (5) is rewritten as 2c (y 2 + c t - tr(WtXt))2 - c(yt - tr(UXt))2 r2((U, Wt) - (U, Wt+1)) Adding the bounds for t = 1, , T, we get 2c L 2 + c MEG(S) - cLU(S) r2((U, W1) - (U, Wt+1)) r2(U, W1), which is equivalent to (6). Assuming LU(S) max and (U, W1) dmax, the bound (6) is tightest when c = r 2dmax/ max. Then we have LMEG(S) - LU(S) r2 maxdmax + r2(U,W 2 1). 3.2 Numerically stable MEG update The MEG update is numerically unstable when the eigenvalues of Wt are around zero. However we can "unwrap" Wt+1 as follows: 1 t Wt+1 = exp(c (^ y ~ tI + log W1 s Z - 2 - ys)Xs), (7) t s=1 where the constant ~ Zt normalizes the trace of Wt+1 to one. As long as the eigen values of W1 are not too small then the computation of log Wt is stable. Note that the update is inde- pendent of the choice of ct R. We incrementally maintain an eigenvalue decomposition of the matrix in the exponent (O(n3) per iteration): t VttVTt = ctI + log W1 - 2 (^ys - ys)Xs), s=1 where the constant ct is chosen so that the maximum eigenvalue of the above is zero. Now Wt+1 = Vt exp(t)VTt /tr(exp(t)). 4 Bregman Projection and DefiniteBoost In this section, we address the following Bregman projection problem2 W = argmin W F (W, W1), tr(W) = 1, tr(WCj ) 0, for j = 1, . . . , n, (8) where the symmetric positive definite matrix W1 of trace one is the initial parameter ma- trix, and C1, . . . , Cn are arbitrary symmetric matrices. Prior knowledge about W is en- coded in the constraints, and the matrix closest to W1 is chosen among the matrices satis- fying all constraints. Tsuda and Noble [13] employed this approach for learning a kernel matrix among graph nodes, and this method can be potentially applied to learn a kernel matrix in other settings (e.g. [14, 11]). The problem (8) is a projection of W1 to the intersection of convex regions defined by the constraints. It is well known that the Bregman projection into the intersection of convex regions can be solved by sequential projections to each region [1]. In the original papers only asymptotic convergence was shown. More recently a connection [4, 7] was made to the AdaBoost algorithm which has an improved convergence analysis [2, 9]. We generalize the latter algorithm and its analysis to symmetric positive definite matrices and call the new algorithm DefiniteBoost. As in the original setting, only approximate projections (Figure 1) are required to show fast convergence. 2Note that if is large then the on-line update (2) becomes a Bregman projection subject to a single equality constraint tr(WXt) = yt. Approximate Figure 1: In (exact) Bregman projections, the intersection Projection of convex sets (i.e., two lines here) is found by iterating pro- jections to each set. We project only approximately, so the projected point does not satisfy the current constraint. Nev- ertheless, global convergence to the optimal solution is guar- anteed via our proofs. Exact Projection Before presenting the algorithm, let us derive the dual problem of (8) by means of Lagrange multipliers , n = argmin log trexp(logW1- jCj), j 0. (9) j=1 See [13] for a detailed derivation of the dual problem. When (8) is feasible, the opti- mal solution is described as W = 1 exp(log W C ) = Z( 1 j ), where Z( ) - nj=1 j tr[exp(log W1 - n C j=1 j j )]. 4.1 Exact Bregman Projections First, let us present the exact Bregman projection algorithm to solve (8). We start from the initial parameter W1. At the t-th step, the most unsatisfied constraint is chosen, jt = argmaxj=1, ,n tr(WtCj). Let us use Ct as the short notation for Cj . Then, the t following Bregman projection with respect to the chosen constraint is solved. Wt+1 = argmin (W, W W t), tr(W) = 1, tr(WCt) 0. (10) By means of a Lagrange multiplier , the dual problem is described as t = argmin tr[exp(log Wt - Ct)], 0. (11) Using the solution of the dual problem, Wt is updated as 1 Wt+1 = exp(log W Z t - tCt) (12) t(t) where the normalization factor is Zt(t) = tr[exp(log Wt -tCt)]. Note that we can use the same numerically stable update as in the previous section. 4.2 Approximate Bregman Projections The solution of (11) cannot be obtained in closed form. However, one can use the following approximate solution: 1 1 + r t/max t t = log , (13) max t - min t 1 + rt/min t when the eigenvalues of Ct lie in the interval [min t , max t ] and rt = tr(WtCt). Since the most unsatisfied constraint is chosen, rt 0 and thus t 0. Although the projection is done only approximately,3 the convergence of the dual objective (9) can be shown using the following upper bound. 3The approximate Bregman projection (with t as in (13) can also be motivated as an online algorithm based on an entropic loss and learning rate one (following Section 3 and [4]). Theorem 4.1 The dual objective (9) is bounded as n T tr explogW1- jCj (rt) (14) j=1 t=1 max min t -t rt max max t -min t rt t -min t where (rt) = 1 - 1 . max - t min t The dual objective is monotonically decreasing, because (rt) 1. Also, since rt corre- sponds to the maximum value among all constraint violations {rj}nj=1, we have (rt) = 1 only if rt = 0. Thus the dual objective continues to decrease until all constraints are satisfied. 4.3 Relation to Boosting When all matrices are diagonal, the DefiniteBoost degenerates to AdaBoost [9]: Let {xi,yi}di=1 be the training samples, where xi Rm and yi {-1,1}. Let h1(x), . . . , hn(x) [-1,1] be the weak hypotheses. For the j-th hypothesis hj(x), let us define Cj = diag(y1hj(x1), . . . , ydhj(xd)). Since |yhj(x)| 1, max/min t = 1 for any t. Setting W1 = I/d, the dual objective (14) is rewritten as d n 1 exp d -yi jhj(xi), i=1 j=1 which is equivalent to the exponential loss function used in AdaBoost. Since Cj and W1 are diagonal, the matrix Wt stays diagonal after the update. If wti = [Wt]ii, the updating formula (12) becomes the AdaBoost update: wt+1,i = wti exp(-tyiht(xi))/Zt(t). The approximate solution of t (13) is described as t = 1 log 1+rt , where r 2 1-r t is the weighted t training error of the t-th hypothesis, i.e. rt = d w i=1 tiyiht(xi). 5 Experiments on Learning Kernels In this section, our technique is applied to learning a kernel matrix from a set of distance measurements. This application is not on-line per se, but it shows nevertheless that the theoretical bounds can be reasonably tight on natural data. When K is a d d kernel matrix among d objects, then the Kij characterizes the similarity between objects i and j. In the feature space, Kij corresponds to the inner product between object i and j, and thus the Euclidean distance can be computed from the entries of the kernel matrix [10]. In some cases, the kernel matrix is not given explicitly, but only a set of distance measurements is available. The data are represented either as (i) quantitative distance values (e.g., the distance between i and j is 0.75), or (ii) qualitative evaluations (e.g., the distance between i and j is small) [14, 13]. Our task is to obtain a positive definite kernel matrix which fits well to the given distance data. On-line kernel learning In the first experiment, we consider the on-line learning scenario in which only one distance example is shown to the learner at each time step. The distance example at time t is described as {at, bt, yt}, which indicates that the squared Euclidean distance between objects at and bt is yt. Let us define a time-developing sequence of kernel matrices as {Wt}Tt=1, and the corresponding points in the feature space as {xti}di=1 (i.e. [Wt]ab = xtaxtb). Then, the total loss incurred by this sequence is T T 2 2 xta = (tr(W t - xtbt - yt tXt) - yt)2, t=1 t=1 1.8 0.45 1.6 0.4 1.4 0.35 1.2 0.3 1 0.25 0.8 Total Loss 0.2 0.6 Classification Error 0.4 0.15 0.2 0.1 0 0.05 0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3 5 Iterations 5 Iterations x 10 x 10 Figure 2: Numerical results of on-line learning. (Left) total loss against the number of iterations. The dashed line shows the loss bound. (Right) classification error of the nearest neighbor classifier using the learned kernel. The dashed line shows the error by the target kernel. where Xt is a symmetric matrix whose (at, at) and (bt, bt) elements are 0.5, (at, bt) and (bt, at) elements are -0.5, and all the other elements are zero. We consider a controlled experiment in which the distance examples are created from a known target kernel matrix. We used a 52 52 kernel matrix among gyrB proteins of bacteria (d = 52). This data contains three bacteria species (see [12] for details). Each distance example is created by randomly choosing one element of the target kernel. The initial parameter was set as W1 = I/d. When the comparison matrix U is set to the target matrix, LU (S) = 0 and max = 0, because all the distance examples are derived from the target matrix. Therefore we choose learning rate = 2, which minimizes the relative loss bound of Lemma 3.2. The total loss of the kernel matrix sequence obtained by the matrix exponential update is shown in Figure 2 (left). In the plot, we have also shown the relative loss bound. The bound seems to give a reasonably tight performance guarantee--it is about twice the actual total loss. To evaluate the learned kernel matrix, the prediction accuracy of bacteria species by the nearest neighbor classifier is calculated (Figure 2, right), where the 52 proteins are randomly divided into 50% training and 50% testing data. The value shown in the plot is the test error averaged over 10 different divisions. It took a large number of iterations ( 2 105) for the error rate to converge to the level of the target kernel. In practice one can often increase the learning rate for faster convergence, but here we chose the small rate suggested by our analysis to check the tightness of the bound. Kernel learning by Bregman projection Next, let us consider a batch learning sce- nario where we have a set of qualitative distance evaluations (i.e. inequality constraints). Given n pairs of similar objects {aj, bj}nj=1, the inequality constraints are constructed as xaj - xbj , j = 1, . . ., n, where is a predetermined constant. If Xj is de- fined as in the previous section and Cj = Xj - I, the inequalities are then rewritten as tr(WCj) 0,j = 1,...,n. The largest and smallest eigenvalues of any Cj are 1 - and -, respectively. As in the previous section, distance examples are generated from the target kernel matrix between gyrB proteins. Setting = 0.2/d, we collected all object pairs whose distance in the feature space is less than to yield 980 inequalities (n = 980). Figure 3 (left) shows the convergence of the dual objective function as proven in Theo- rem 4.1. The convergence was much faster than the previous experiment, because, in the batch setting, one can choose the most unsatisfied constraint, and optimize the step size as well. Figure 3 (right) shows the classification error of the nearest neighbor classifier. As opposed to the previous experiment, the error rate is higher than that of the target kernel matrix, because substantial amount of information is lost by the conversion to inequality constraints. 55 0.8 50 0.7 45 0.6 40 0.5 35 0.4 Dual Obj 30 0.3 Classification Error 25 0.2 20 0.1 15 0 0 50 100 150 200 250 300 0 50 100 150 200 250 300 Iterations Iterations Figure 3: Numerical results of Bregman projection. (Left) convergence of the dual objective function. (Right) classification error of the nearest neighbor classifier using the learned kernel. 6 Conclusion We motivated and analyzed a new update for symmetric positive matrices using the von Neumann divergence. We showed that the standard bounds for on-line learning and Boost- ing generalize to the case when the parameters are a symmetric positive definite matrix (of trace one) instead of a probability vector. As in quantum physics, the eigenvalues act as probabilities. Acknowledgment We would like to thank B. Sch olkopf, M. Kawanabe, J. Liao and W.S. Noble for fruitful discussions. M.W. was supported by NSF grant CCR 9821087 and UC Discovery grant LSIT02-10110. K.T. and G.R. gratefully acknowledge partial support from the PASCAL Network of Excellence (EU #506778). Part of this work was done while all three authors were visiting the National ICT Australia in Canberra.

This letter analyzes the Fisher kernel from a statistical point of view. The Fisher kernel is a particularly interesting method for constructing a model of the posterior probability that makes intelligent use of unlabeled data (i.e., of the underlying data density). It is important to analyze and ultimately understand the statistical properties of the Fisher kernel. To this end, we first establish sufficient conditions that the constructed posterior model is realizable (i.e., it contains the true distribution). Realizability immediately leads to consistency results. Subsequently, we focus on an asymptotic analysis of the generalization error, which elucidates the learning curves of the Fisher kernel and how unlabeled data contribute to learning. We also point out that the squared or log loss is theoretically preferable-because both yield consistent estimators-to other losses such as the exponential loss, when a linear classifier is used together with the Fisher kernel. Therefore, this letter underlines that the Fisher kernel should be viewed not as a heuristics but as a powerful statistical tool with well-controlled statistical properties.

A common way of image denoising is to project a noisy image to the sub- space of admissible images made for instance by PCA. However, a major drawback of this method is that all pixels are updated by the projection, even when only a few pixels are corrupted by noise or occlusion. We pro- pose a new method to identify the noisy pixels by (cid:1) 1-norm penalization and update the identified pixels only. The identification and updating of noisy pixels are formulated as one linear program which can be solved efficiently. Especially, one can apply the ν-trick to directly specify the fraction of pixels to be reconstructed. Moreover, we extend the linear program to be able to exploit prior knowledge that occlusions often ap- pear in contiguous blocks (e.g. sunglasses on faces). The basic idea is to penalize boundary points and interior points of the occluded area dif- ferently. We are able to show the ν-property also for this extended LP leading a method which is easy to use. Experimental results impressively demonstrate the power of our approach.

MOTIVATION: Kernel methods such as support vector machines require a kernel function between objects to be defined a priori. Several works have been done to derive kernels from probability distributions, e.g., the Fisher kernel. However, a general methodology to design a kernel is not fully developed. RESULTS: We propose a reasonable way of designing a kernel when objects are generated from latent variable models (e.g., HMM). First of all, a joint kernel is designed for complete data which include both visible and hidden variables. Then a marginalized kernel for visible data is obtained by taking the expectation with respect to hidden variables. We will show that the Fisher kernel is a special case of marginalized kernels, which gives another viewpoint to the Fisher kernel theory. Although our approach can be applied to any object, we particularly derive several marginalized kernels useful for biological sequences (e.g., DNA and proteins). The effectiveness of marginalized kernels is illustrated in the task of classifying bacterial gyrase subunit B (gyrB) amino acid sequences.

Recently the Fisher score (or the Fisher kernel) is increasingly used as a feature extractor for classification problems. The Fisher score is a vector of parameter derivatives of loglikelihood of a probabilistic model. This paper gives a theoretical analysis about how class information is pre- served in the space of the Fisher score, which turns out that the Fisher score consists of a few important dimensions with class information and many nuisance dimensions. When we perform clustering with the Fisher score, K-Means type methods are obviously inappropriate because they make use of all dimensions. So we will develop a novel but simple clus- tering algorithm specialized for the Fisher score, which can exploit im- portant dimensions. This algorithm is successfully tested in experiments with artificial data and real data (amino acid sequences). as follows: Let us assume that a probabilistic model   parameter estimate 

Recently, Jaakkola and Haussler (1999) proposed a method for constructing kernel functions from probabilistic models. Their so-called Fisher kernel has been combined with discriminative classifiers such as support vector machines and applied successfully in, for example, DNA and protein analysis. Whereas the Fisher kernel is calculated from the marginal log-likelihood, we propose the TOP kernel derived; from tangent vectors of posterior log-odds. Furthermore, we develop a theoretical framework on feature extractors from probabilistic models and use it for analyzing the TOP kernel. In experiments, our new discriminative TOP kernel compares favorably to the Fisher kernel.

Nonquadratic regularizers, in particular the l(1) norm regularizer can yield sparse solutions that generalize well. In this work we propose the generalized subspace information criterion (GSIC) that allows to predict the generalization error for this useful family of regularizers. We show that under some technical assumptions GSIC is an asymptotically unbiased estimator of the generalization error. GSIC is demonstrated to have a good performance in experiments with the l(1) norm regularizer as we compare with the network information criterion (NIC) and cross- validation in relatively large sample cases. However in the small sample case, GSIC tends to fail to capture the optimal model due to its large variance. Therefore, also a biased version of GSIC is introduced,which achieves reliable model selection in the relevant and challenging scenario of high-dimensional data and few samples.

Machine recognition of hand-printed Japanese characters has been an area of great interest for many years. A major problem of this classification task is the huge number of different characters. Applying standard "state-of-the-art" techniques, such as SVM, to multi-class problems of this kind imposes severe problems of both a conceptual and technical nature: (i) separating one class from all others may be an unnecessarily hard problem; and (ii) solving these subproblems can impose unacceptably high computational costs. In this paper, a new approach to Japanese character recognition is presented that successfully overcomes these shortcomings. It is based on a pairwise coupling procedure for probabilistic two-class kernel classifiers. Experimental results for Hiragana recognition effectively demonstrate that our method attains an excellent level of prediction accuracy while imposing very low computational costs.

Recently, Jaakkola and Haussler proposed a method for construct(cid:173) ing kernel functions from probabilistic models. Their so called "Fisher kernel" has been combined with discriminative classifiers such as SVM and applied successfully in e.g. DNA and protein analysis. Whereas the Fisher kernel (FK) is calculated from the marginal log-likelihood, we propose the TOP kernel derived from Tangent vectors Of Posterior log-odds. Furthermore we develop a theoretical framework on feature extractors from probabilistic models and use it for analyzing FK and TOP. In experiments our new discriminative TOP kernel compares favorably to the Fisher kernel.

This paper provides an introduction to support vector machines, kernel Fisher discriminant analysis, and kernel principal component analysis, as examples for successful kernel-based learning methods. We first give a short background about Vapnik-Chervonenkis theory and kernel feature spaces and then proceed to kernel based learning in supervised and unsupervised scenarios including practical and algorithmic considerations. We illustrate the usefulness of kernel algorithms by discussing applications such as optical character recognition and DNA analysis.

This is an elementary research into assigning color values to voxels of multi-channel magnetic resonance imaging (MRI) volume data. The MRI volume data sets obtained under different scanning conditions are transformed into components by independent component analysis (ICA), which enhances the physical characteristics of the tissue. The transfer functions for generating color values from the independent components are obtained by using a radial basis function network, a kind of neural net, by training the network with sample data chosen from the Visible Human female data set (VHF). The resultant color volume data sets correspond well with the full-color cross-sections of the Visible Human data sets.

In this paper, we construct a learning method of stochastic perceptron based on semiparametric inference, and show that this method produces large margin solutions. In semiparametric inference, the parameters are divided into structural parameters which are to be estimated and nuisance parameters in which we do not have any interest. Here, the weight vector of perceptron is defined as structural parameters and the steepness of transfer function is defined as a nuisance parameter. Usually, rough estimate is substituted to nuisance parameters and only structural parameters are estimated. To compensate the estimation error caused by rough estimate, an additional term is added to the derivative of likelihood. We will show that this additional term is related to the regularization term which causes large margin solutions. This work suggests that the success of large margin classifiers can be attributed to semiparametric inference.

To improve the performance of subspace classifier, it is effective to reduce the dimensionality of the intersections between subspaces. For this purpose, the feature space is mapped implicitly to a high dimensional reproducing kernel Hilbert space and the subspace classifier is applied in this space. As a result of Hiragana recognition experiment, our classifier outperformed the conventional subspace classifier.

When the dimensionality of the feature space increases and takes beyond a certain point, the classification performance of a parametric classifier begins to deteriorate. This is because the number of parameters of the classifier depends on the dimensionality and gets too large in a high dimensional space. To obtain the density value at the point of interest without deriving the parameters, we propose a nonparametric Gaussian density model, where the sum of the log-density values at two points is described, without any parameter, by a function of the distance between the two points. The density value at the point of interest is estimated from the distance between the point and each of the training sample points using this model. We will empirically show that the density values estimated by the nonparametric model are more accurate than those by the parametric model in a high dimensional space. In a character recognition experiment, our nonparametric classifier achieved higher classification accuracy in comparison with the parametric classifier.

Document clustering is a powerful tool for browsing throughout a document database. Similar documents are gathered into several clusters and a representative document of each cluster is shown to users. To make users infer the content of the database from several representatives, the documents must be separated into tight clusters, in which documents are connected with high similarities. At the same time, clustering must be fast for user interaction. We propose an O(n/sup 2/) time, O(n) space cluster extraction method. It is faster than the ordinal clustering methods, and its clusters compare favorably with those produced by Complete Link for tightness. When we deal with OCR-ed documents, term loss caused by recognition faults can change similarities between documents. We also examined the effect of recognition faults to the performance of document clustering.