Solving the model counting problem #SAT, asking for the number of satisfying assignments of a propositional formula, has been explored intensively and has gathered its own community. While most existing solvers are based on knowledge compilation, another promising approach is through contraction in tensor hypernetworks. We perform a theoretical proof-complexity analysis of this approach. For this, we design two new tensor-based proof systems that we show to tightly correspond to tensor-based #SAT solving. We determine the simulation order of #SAT proof systems and prove exponential separations between the systems. This sheds light on the relative performance of different #SAT solving approaches.

Solutions of symbolic regression problems are expressions that are composed of input variables and operators from a finite set of function symbols. One measure for evaluating symbolic regression algorithms is their ability to recover formulae, up to symbolic equivalence, from finite samples. Not unexpectedly, the recovery problem becomes harder when the formula gets more complex, that is, when the number of variables and operators gets larger. Variables in naturally occurring symbolic formulas often appear only in fixed combinations. This can be exploited in symbolic regression by substituting one new variable for the combination, effectively reducing the number of variables. However, finding valid substitutions is challenging. Here, we address this challenge by searching over the expression space of small substitutions and testing for validity. The validity test is reduced to a test of functional dependence. The resulting iterative dimension reduction procedure can be used with any symbolic regression approach. We show that it reliably identifies valid substitutions and significantly boosts the performance of different types of state-of-the-art symbolic regression algorithms.

Solving systems of ordinary differential equations (ODEs) is essential when it comes to understanding the behavior of dynamical systems. Yet, automated solving remains challenging, in particular for nonlinear systems. Computer algebra systems (CASs) provide support for solving ODEs by first simplifying them, in particular through the use of Lie point symmetries. Finding these symmetries is, however, itself a difficult problem for CASs. Recent works in symbolic regression have shown promising results for recovering symbolic expressions from data. Here, we adapt search-based symbolic regression to the task of finding generators of Lie point symmetries. With this approach, we can find symmetries of ODEs that existing CASs cannot find.

Einsum is a declarative language for tensor expressions that specifies an output tensor in terms of several input tensors. However, it does not specify how to compute the output tensor from the input tensors. A typical computational backend for the einsum language comprises two parts: First, a contraction path algorithm that breaks down an einsum expression into a sequence of binary tensor contractions. Second, the execution of the binary contractions. For efficient binary contractions, the data layout of the tensors must be optimized. So far, the computation of contraction paths and the optimization of the data layout for single, that is, local, binary tensor contractions have been studied in isolation. For optimizing the overall execution times of einsum expressions, we introduce Einsum Tree IR, an intermediate representation for globally optimizing the data layout for a given contraction path. We illustrate the effectiveness of the approach on a state-of-the-art Arm server processor, an x86 server processor, and an x86 desktop system.

Einsum expressions specify an output tensor in terms of several input tensors. They offer a simple yet expressive abstraction for many computational tasks in artificial intelligence and beyond. However, evaluating einsum expressions poses hard algorithmic problems that depend on the representation of the tensors. Two popular representations are multidimensional arrays and coordinate lists. The latter is a more compact representation for sparse tensors, that is, tensors where a significant proportion of the entries are zero. So far, however, most of the popular einsum implementations use the multidimensional array representation for tensors. Here, we show on a non-trivial example that, when evaluating einsum expressions, coordinate lists can be exponentially more efficient than multidimensional arrays. In practice, however, coordinate lists can also be significantly less efficient than multidimensional arrays, but it is hard to decide from the input tensors whether this will be the case. Sparsity evolves dynamically in intermediate tensors during the evaluation of an einsum expression. Therefore, we introduce a hybrid solution where the representation is switched on the fly from multidimensional arrays to coordinate lists depending on the sparsity of the remaining tensors. In our experiments on established benchmark einsum expressions, the hybrid solution is consistently competitive with or outperforms the better of the two static representations.

Many decision and optimization problems have natural extensions as counting problems. The best known example is the Boolean satisfiability problem (SAT), where we want to count the satisfying assignments of truth values to the variables, which is known as the #SAT problem. Likewise, for discrete optimization problems, we want to count the states on which the objective function attains the optimal value. Both SAT and discrete optimization can be formulated as selective marginalize a product function (MPF) queries. Here, we show how general selective MPF queries can be extended for model counting. MPF queries are encoded as tensor hypernetworks over suitable semirings that can be solved by generic tensor hypernetwork contraction algorithms. Our model counting extension is again an MPF query, on an extended semiring, that can be solved by the same contraction algorithms. Model counting is required for uniform model sampling. We show how the counting extension can be further extended for model sampling by constructing yet another semiring. We have implemented the model counting and sampling extensions. Experiments show that our generic approach is competitive with the state of the art in model counting and model sampling.

Modern artificial intelligence and machine learning workflows rely on efficient tensor libraries. However, tuning tensor libraries without considering the actual problems they are meant to execute can lead to a mismatch between expected performance and the actual performance. Einsum libraries are tuned to efficiently execute tensor expressions with only a few, relatively large, dense, floating-point tensors. But, practical applications of einsum cover a much broader range of tensor expressions than those that can currently be executed efficiently. For this reason, we have created a benchmark dataset that encompasses this broad range of tensor expressions, allowing future implementations of einsum to build upon and be evaluated against. In addition, we also provide generators for einsum expressions and converters to einsum expressions in our repository, so that additional data can be generated as needed. The benchmark dataset, the generators and converters are released openly and are publicly available at https://benchmark.einsum.org.

Gaussian mixture models are classical but still popular machine learning models. An appealing feature of Gaussian mixture models is their tractability, that is, they can be learned efficiently and exactly from data, and also support efficient exact inference queries like soft clustering data points. Only seemingly simple, Gaussian mixture models can be hard to understand. There are at least four aspects to understanding Gaussian mixture models, namely, understanding the whole distribution, its individual parts (mixture components), the relationships between the parts, and the interplay of the whole and its parts. In a structured literature review of applications of Gaussian mixture models, we found the need for supporting all four aspects. To identify candidate visualizations that effectively aid the user needs, we structure the available design space along three different representations of Gaussian mixture models, namely as functions, sets of parameters, and sampling processes. From the design space, we implemented three design concepts that visualize the overall distribution together with its components. Finally, we assessed the practical usefulness of the design concepts with respect to the different user needs in expert interviews and an insight-based user study.

Capsule neural networks replace simple, scalar-valued neurons with vector-valued capsules. They are motivated by the pattern recognition system in the human brain, where complex objects are decomposed into a hierarchy of simpler object parts. Such a hierarchy is referred to as a parse-tree. Conceptually, capsule neural networks have been defined to mimic this behavior. The capsule neural network (CapsNet), by Sabour, Frosst, and Hinton, is the first actual implementation of the conceptual idea of capsule neural networks. CapsNets achieved state-of-the-art performance on simple image recognition tasks with fewer parameters and greater robustness to affine transformations than comparable approaches. This sparked extensive follow-up research. However, despite major efforts, no work was able to scale the CapsNet architecture to more reasonable-sized datasets. Here, we provide a reason for this failure and argue that it is most likely not possible to scale CapsNets beyond toy examples. In particular, we show that the concept of a parse-tree, the main idea behind capsule neuronal networks, is not present in CapsNets. We also show theoretically and experimentally that CapsNets suffer from a vanishing gradient problem that results in the starvation of many capsules during training.

Computational problems ranging from artificial intelligence to physics require efficient computations of large tensor expressions. These tensor expressions can often be represented in Einstein notation. To evaluate tensor expressions in Einstein notation, that is, for the actual Einstein summation, usually external libraries are used. Surprisingly, Einstein summation operations on tensors fit well with fundamental SQL constructs. We show that by applying only four mapping rules and a simple decomposition scheme using common table expressions, large tensor expressions in Einstein notation can be translated to portable and efficient SQL code. The ability to execute large Einstein summation queries opens up new possibilities to process data within SQL. We demonstrate the power of Einstein summation queries on four use cases, namely querying triplestore data, solving Boolean satisfiability problems, performing inference in graphical models, and simulating quantum circuits. The performance of Einstein summation queries, however, depends on the query engine implemented in the database system. Therefore, supporting efficient Einstein summation computations in database systems presents new research challenges for the design and implementation of query engines.

The Gaussian mixture model (GMM) describes the distribution of random variables from several different populations. GMMs have widespread applications in probability theory, statistics, machine learning for unsupervised cluster analysis and topic modeling, as well as in deep learning pipelines. So far, few efforts have been made to explore the underlying point distribution in combination with the GMMs, in particular when the data becomes high-dimensional and when the GMMs are composed of many Gaussians. We present an analysis tool comprising various GPU-based visualization techniques to explore such complex GMMs. To facilitate the exploration of high-dimensional data, we provide a novel navigation system to analyze the underlying data. Instead of projecting the data to 2D, we utilize interactive 3D views to better support users in understanding the spatial arrangements of the Gaussian distributions. The interactive system is composed of two parts: (1) raycasting-based views that visualize cluster memberships, spatial arrangements, and support the discovery of new modes. (2) overview visualizations that enable the comparison of Gaussians with each other, as well as small multiples of different choices of basis vectors. Users are supported in their exploration with customization tools and smooth camera navigations. Our tool was developed and assessed by five domain experts, and its usefulness was evaluated with 23 participants. To demonstrate the effectiveness, we identify interesting features in several data sets.

Probabilistic programming is a powerful means for formally specifying machine learning models. The inference engine of a probabilistic programming environment can be used for serving complex queries on these models. Most of the current research in probabilistic programming is dedicated to the design and implementation of highly efficient inference engines. Much less research aims at making the power of these inference engines accessible to non-expert users. Probabilistic programming means writing code. Yet many potential users from promising application areas such as the social sciences lack programming skills. This prompted recent efforts in synthesizing probabilistic programs directly from data. However, working with synthesized programs still requires the user to read, understand, and write some code, for instance, when invoking the inference engine for answering queries. Here, we present an interactive visual approach to synthesizing and querying probabilistic programs that does not require the user to read or write code.

Constrained optimization problems arise frequently in classical machine learning. There exist frameworks addressing constrained optimization, for instance, CVXPY and GENO. However, in contrast to deep learning frameworks, GPU support is limited. Here, we extend the GENO framework to also solve constrained optimization problems on the GPU. The framework allows the user to specify constrained optimization problems in an easy-to-read modeling language. A solver is then automatically generated from this specification. When run on the GPU, the solver outperforms state-of-the-art approaches like CVXPY combined with a GPU-accelerated solver such as cuOSQP or SCS by a few orders of magnitude.

Deep learning models for different NLP tasks often rely on pre-trained word embeddings, that is, vector representations of words. Therefore, it is crucial to evaluate pre-trained word embeddings independently of downstream tasks. Such evaluations try to assess whether the geometry induced by a word embedding captures connections made in natural language, such as, analogies, clustering of words, or word similarities. Here, traditionally, similarity is measured by comparison to human judgment. However, explicitly annotating word pairs with similarity scores by surveying humans is expensive. We tackle this problem by formulating a similarity measure that is based on an agent for routing the Wikipedia hyperlink graph. In this graph, word similarities are implicitly encoded by edges between articles. We show on the English Wikipedia that our measure correlates well with a large group of traditional similarity measures, while covering a much larger proportion of words and avoiding explicit human labeling. Moreover, since Wikipedia is available in more than 300 languages, our measure can easily be adapted to other languages, in contrast to traditional similarity measures.

The Hessian of a differentiable convex function is positive semidefinite. Therefore, checking the Hessian of a given function is a natural approach to certify convexity. However, implementing this approach is not straightforward, since it requires a representation of the Hessian that allows its analysis. Here, we implement this approach for a class of functions that is rich enough to support classical machine learning. For this class of functions, it was recently shown how to compute computational graphs of their Hessians. We show how to check these graphs for positive-semidefiniteness. We compare our implementation of the Hessian approach with the well-established disciplined convex programming (DCP) approach and prove that the Hessian approach is at least as powerful as the DCP approach for differentiable functions. Furthermore, we show for a state-of-the-art implementation of the DCP approach that the Hessian approach is actually more powerful, that is, it can certify the convexity of a larger class of differentiable functions.

Abstract Recent works showed that implementations of quicksort using vector CPU instructions can outperform the non‐vectorized algorithms in widespread use. However, these implementations are typically single‐threaded, implemented for a particular instruction set, and restricted to a small set of key types. We lift these three restrictions: our proposedvqsortalgorithm integrates into the state‐of‐the‐art parallel sorter , with a geometric mean speedup of 1.59. The same implementation works on seven instruction sets (including SVE and RISC‐V V) across four platforms. It also supports floating‐point and 16–128 bit integer keys. To the best of our knowledge, this is the fastest sort for large arrays of non‐tuple keys on CPUs, up to 20 times as fast as the sorting algorithms implemented in standard libraries. This article focuses on the practical engineering aspects enabling the speed and portability, which we have not yet seen demonstrated for a quicksort implementation. Furthermore, we introduce compact and transpose‐free sorting networks for in‐register sorting of small arrays, and a vector‐friendly pivot sampling strategy that is robust against adversarial input.

By skeptics and undecided we refer to nodes in clustered social networks that cannot be assigned easily to any of the clusters. Such nodes are typically found either at the interface between clusters (the undecided) or at their boundaries (the skeptics). Identifying these nodes is relevant in marketing applications like voter targeting, because the persons represented by such nodes are often more likely to be affected in marketing campaigns than nodes deeply within clusters. So far this identification task is not as well studied as other network analysis tasks like clustering, identifying central nodes, and detecting motifs. We approach this task by deriving novel geometric features from the network structure that naturally lend themselves to an interactive visual approach for identifying interface and boundary nodes.

Topic models are characterized by a latent class variable that represents the different topics. Traditionally, their observable variables are modeled as discrete variables like, for instance, in the prototypical latent Dirichlet allocation (LDA) topic model. In LDA, words in text documents are encoded by discrete count vectors with respect to some dictionary. The classical approach for learning topic models optimizes a likelihood function that is non-concave due to the presence of the latent variable. Hence, this approach mostly boils down to using search heuristics like the EM algorithm for parameter estimation. Recently, it was shown that topic models can be learned with strong algorithmic and statistical guarantees through Pearson's method of moments. Here, we extend this line of work to topic models that feature discrete as well as continuous observable variables (features). Moving beyond discrete variables as in LDA allows for more sophisticated features and a natural extension of topic models to other modalities than text, like, for instance, images. We provide algorithmic and statistical guarantees for the method of moments applied to the extended topic model that we corroborate experimentally on synthetic data. We also demonstrate the applicability of our model on real-world document data with embedded images that we preprocess into continuous state-of-the-art feature vectors.

It has long been known that principal component analysis (PCA) is not robust with respect to gross data corruption. This has been addressed by robust principal component analysis (RPCA). The first computationally tractable definition of RPCA decomposes a data matrix into a low-rank and a sparse component. The low-rank component represents the principal components, while the sparse component accounts for the data corruption. Previous works consider the corruption of individual entries or whole columns of the data matrix. In contrast, we consider a more general form of data corruption that affects groups of measurements. We show that the decomposition approach remains computationally tractable and allows the exact recovery of the decomposition when only the corrupted data matrix is given. Experiments on synthetic data corroborate our theoretical findings, and experiments on several real-world datasets from different domains demonstrate the wide applicability of our generalized approach.

Modern CPUs provide single instruction-multiple data (SIMD) instructions. SIMD instructions process several elements of a primitive data type simultaneously in fixed-size vectors. Classical sorting algorithms are not directly expressible in SIMD instructions. Accelerating sorting algorithms with SIMD instruction is therefore a creative endeavor. A promising approach for sorting with SIMD instructions is to use sorting networks for small arrays and Quicksort for large arrays. In this paper we improve vectorization techniques for sorting networks and Quicksort. In particular, we show how to use the full capacity of vector registers in sorting networks and how to make vectorized Quicksort robust with respect to different key distributions. To demonstrate the performance of our techniques we implement an in-place hybrid sorting algorithm for the data type int with AVX2 intrinsics. Our implementation is at least 30% faster than state-of-the-art high-performance sorting alternatives.

Most problems from classical machine learning can be cast as an optimization problem. We introduce GENO (GENeric Optimization), a framework that lets the user specify a constrained or unconstrained optimization problem in an easy-to-read modeling language. GENO then generates a solver, i.e., Python code, that can solve this class of optimization problems. The generated solver is usually as fast as hand-written, problem-specific, and well-engineered solvers. Often the solvers generated by GENO are faster by a large margin compared to recently developed solvers that are tailored to a specific problem class.An online interface to our framework can be found at http://www.geno-project.org.

Computing derivatives of tensor expressions, also known as tensor calculus, is a fundamental task in machine learning. A key concern is the efficiency of evaluating the expressions and their derivatives that hinges on the representation of these expressions. Recently, an algorithm for computing higher order derivatives of tensor expressions like Jacobians or Hessians has been introduced that is a few orders of magnitude faster than previous state-of-the-art approaches. Unfortunately, the approach is based on Ricci notation and hence cannot be incorporated into automatic differentiation frameworks like TensorFlow, PyTorch, autograd, or JAX that use the simpler Einstein notation. This leaves two options, to either change the underlying tensor representation in these frameworks or to develop a new, provably correct algorithm based on Einstein notation. Obviously, the first option is impractical. Hence, we pursue the second option. Here, we show that using Ricci notation is not necessary for an efficient tensor calculus and develop an equally efficient method for the simpler Einstein notation. It turns out that turning to Einstein notation enables further improvements that lead to even better efficiency.

Measurement is at the core of scientific discovery. However, some quantities, such as economic behavior or intelligence, do not allow for direct measurement. They represent latent constructs that require surrogate measurements. In other scenarios, non-observed quantities can influence the variables of interest. In either case, models with latent variables are needed. Here, we investigate fused latent and graphical models that exhibit continuous latent variables and discrete observed variables. These models are characterized by a decomposition of the pairwise interaction parameter matrix into a group-sparse component of direct interactions and a low-rank component of indirect interactions due to the latent variables. We first investigate when such a decomposition is identifiable. Then, we show that fused latent and graphical models can be recovered consistently from data in the high-dimensional setting. We support our theoretical findings with experiments on synthetic and real-world data from polytomous item response theory studies.

Computing derivatives of tensor expressions, also known as tensor calculus, is a fundamental task in machine learning. A key concern is the efficiency of evaluating the expressions and their derivatives that hinges on the representation of these expressions. Recently, an algorithm for computing higher order derivatives of tensor expressions like Jacobians or Hessians has been introduced that is a few orders of magnitude faster than previous state-of-the-art approaches. Unfortunately, the approach is based on Ricci notation and hence cannot be incorporated into automatic differentiation frameworks from deep learning like TensorFlow, PyTorch, autograd, or JAX that use the simpler Einstein notation. This leaves two options, to either change the underlying tensor representation in these frameworks or to develop a new, provably correct algorithm based on Einstein notation. Obviously, the first option is impractical. Hence, we pursue the second option. Here, we show that using Ricci notation is not necessary for an efficient tensor calculus and develop an equally efficient method for the simpler Einstein notation. It turns out that turning to Einstein notation enables further improvements that lead to even better efficiency. The methods that are described in this paper for computing derivatives of matrix and tensor expressions have been implemented in the online tool www.MatrixCalculus.org .

Regularized loss minimization, where a statistical model is obtained from minimizing the sum of a loss function and weighted regularization terms, is still in widespread use in machine learning. The statistical performance of the resulting models depends on the choice of weights (regularization parameters) that are typically tuned by cross-validation. For finding the best regularization parameters, the regularized minimization problem needs to be solved for the whole parameter domain. A practically more feasible approach is covering the parameter domain with approximate solutions of the loss minimization problem for some prescribed approximation accuracy. The problem of computing such a covering is known as the approximate solution gamut problem. Existing algorithms for the solution gamut problem suffer from several problems. For instance, they require a grid on the parameter domain whose spacing is difficult to determine in practice, and they are not generic in the sense that they rely on problem specific plug-in functions. Here, we show that a well-known algorithm from vector optimization, namely the Benson algorithm, can be used directly for computing approximate solution gamuts while avoiding the problems of existing algorithms. Experiments for the Elastic Net on real world data sets demonstrate the effectiveness of Benson’s algorithm for regularization parameter tracking.

Ising models describe the joint probability distribution of a vector of binary feature variables. Typically, not all the variables interact with each other and one is interested in learning the presumably sparse network structure of the interacting variables. However, in the presence of latent variables, the conventional method of learning a sparse model might fail. This is because the latent variables induce indirect interactions of the observed variables. In the case of only a few latent conditional {Gaussian} variables these spurious interactions contribute an additional low-rank component to the interaction parameters of the observed Ising model. Therefore, we propose to learn a sparse + low-rank decomposition of the parameters of an {Ising} model using a convex regularized likelihood problem. We show that the same problem can be obtained as the dual of a maximum-entropy problem with a new type of relaxation, where the sample means collectively need to match the expected values only up to a given tolerance. The solution to the convex optimization problem has consistency properties in the high-dimensional setting, where the number of observed binary variables and the number of latent conditional {Gaussian} variables are allowed to grow with the number of training samples.

Many machine learning methods entail minimizing a loss-function that is the sum of the losses for each data point. The form of the loss function is exploited algorithmically, for instance in stochastic gradient descent (SGD) and in the alternating direction method of multipliers (ADMM). However, there are also machine learning methods where the entailed optimization problem features the data points not in the objective function but in the form of constraints, typically one constraint per data point. Here, we address the problem of solving convex optimization problems with many convex constraints. Our approach is an extension of ADMM. The straightforward implementation of ADMM for solving constrained optimization problems in a distributed fashion solves constrained subproblems on different compute nodes that are aggregated until a consensus solution is reached. Hence, the straightforward approach has three nested loops: one for reaching consensus, one for the constraints, and one for the unconstrained problems. Here, we show that solving the costly constrained subproblems can be avoided. In our approach, we combine the ability of ADMM to solve convex optimization problems in a distributed setting with the ability of the augmented Lagrangian method to solve constrained optimization problems. Consequently, our algorithm only needs two nested loops. We prove that it inherits the convergence guarantees of both ADMM and the augmented Lagrangian method. Experimental results corroborate our theoretical findings.

Latent variable graphical models are an extension of Gaussian graphical models that decompose the precision matrix into a sparse and a low-rank component. These models can be learned with theoretical guarantees from data via a semidefinite program. This program features two regularization terms, one for promoting sparsity and one for promoting a low rank. In practice, however, it is not straightforward to learn a good model since the model highly depends on the regularization parameters that control the relative weight of the loss function and the two regularization terms. Selecting good regularization parameters can be modeled as a bi-level optimization problem, where the upper level optimizes some form of generalization error and the lower level provides a description of the solution gamut. The solution gamut is the set of feasible solutions for all possible values of the regularization parameters. In practice, it is often not feasible to describe the solution gamut efficiently. Hence, algorithmic schemes for approximating solution gamuts have been devised. One such scheme is Benson's generic vector optimization algorithm that comes with approximation guarantees. So far Benson's algorithm has not been used in conjunction with semidefinite programs like the latent variable graphical Lasso. Here, we develop an adaptive variant of Benson's algorithm for the semidefinite case and show that it keeps the known approximation and run time guarantees. Furthermore, Benson's algorithm turns out to be practically more efficient for the latent variable graphical model than the existing solution gamut approximation scheme on a wide range of data sets.

Although optimization is the longstanding, algorithmic backbone of machine learning new models still require the time-consuming implementation of new solvers. As a result, there are thousands of implementations of optimization algorithms for machine learning problems. A natural question is, if it is always necessary to implement a new solver, or is there one algorithm that is sufficient for most models. Common belief suggests that such a one-algorithm-fits-all approach cannot work, because this algorithm cannot exploit model specific structure. At least, a generic algorithm cannot be efficient and robust on a wide variety of problems. Here, we challenge this common belief. We have designed and implemented the optimization framework GENO (GENeric Optimization) that combines a modeling language with a generic solver. GENO takes the declaration of an optimization problem and generates a solver for the specified problem class. The framework is flexible enough to encompass most of the classical machine learning problems. We show on a wide variety of classical but also some recently suggested problems that the automatically generated solvers are (1) as efficient as well engineered, specialized solvers, (2) more efficient by a decent margin than recent state-of-the-art solvers, and (3) orders of magnitude more efficient than classical modeling language plus solver approaches.

Abstract The development of custom interactive visualization tools for specific domains and applications has been made much simpler recently by a surge of visualization tools, libraries and frameworks. Most of these tools are developed for classical data science applications, where a user is supported in analyzing measured or simulated data. But recently, there has also been an increasing interest in visual support for understanding machine learning algorithms and frameworks, especially for deep learning. Many, if not most, of the visualization support for (deep) learning addresses the developer of the learning system and not the end user (data scientist). Here we show on a specific example, namely the development of a matrix calculus algorithm, that supporting visualizations can also greatly benefit the development of algorithms in classical domains like in our case computer algebra. The idea is similar to visually supporting the understanding of learning algorithms, namely provide the developer with an interactive, visual tool that provides insights into the workings and, importantly, also into the failures of the algorithm under development. Developing visualization support for matrix calculus development went similar as the development of more traditional visual support systems for data analysts. First, we had to acquaint ourselves with the problem, its language and challenges by talking to the core developer of the matrix calculus algorithm. Once we understood the challenge, it was fairly easy to develop visual support that streamlined the development of the matrix calculus algorithm significantly.

Optimization is an integral part of most machine learning systems and most numerical optimization schemes rely on the computation of derivatives. Therefore, frameworks for computing derivatives are an active area of machine learning research. Surprisingly, as of yet, no existing framework is capable of computing higher order matrix and tensor derivatives directly. Here, we close this fundamental gap and present an algorithmic framework for computing matrix and tensor derivatives that extends seamlessly to higher order derivatives. The framework can be used for symbolic as well as for forward and reverse mode automatic differentiation. Experiments show a speedup between one and four orders of magnitude over state-of-the-art frameworks when evaluating higher order derivatives.

Abstract Scatter plots are mostly used for correlation analysis, but are also a useful tool for understanding the distribution of high‐dimensional point cloud data. An important characteristic of such distributions are clusters, and scatter plots have been used successfully to identify clusters in data. Another characteristic of point cloud data that has received less attention so far are regions that contain no or only very few data points. We show that augmenting scatter plots by projections of flow lines along the gradient vector field of the distance function to the point cloud reveals such empty regions or voids. The augmented scatter plots, that we call sclow plots, enable a much better understanding of the geometry underlying the point cloud than traditional scatter plots, and by that support tasks like dimension inference, detecting outliers, or identifying data points at the interface between clusters. We demonstrate the feasibility of our approach on synthetic and real world data sets.

Many machine learning methods are given as parameterized optimization problems. Important examples of such parameters are regularization- and kernel hyperparameters. These parameters have to be tuned carefully since the choice of their values can have a significant impact on the statistical performance of the learning methods. In most cases the parameter space does not carry much structure and parameter tuning essentially boils down to exploring the whole parameter space. The case when there is only one parameter received quite some attention over the years. First, algorithms for tracking an optimal solution for several machine learning optimization problems over regularization- and hyperparameter intervals had been developed, but since these algorithms can suffer from numerical problems more robust and efficient approximate path tracking algorithms have been devised and analyzed recently. By now approximate path tracking algorithms are known for regularization-and kernel hyperparameter paths with optimal path complexities that depend only on the prescribed approximation error. Here we extend the work on approximate path tracking algorithms with approximation guarantees to multi-dimensional parameter domains. We show a lower bound on the complexity of approximately exploring a multi-dimensional parameter domain that is the product of the corresponding path complexities. We also show a matching upper bound that can be turned into a theoretically and practically efficient algorithm. Experimental results for kernelized support vector machines and the elastic net confirm the theoretical complexity analysis.

We want to sketch the support of a probability measure on Euclidean space from samples that have been drawn from the measure. This problem is closely related to certain manifold learning problems, where one assumes that the sample points are drawn from a manifold that is embedded in Euclidean space. Here we propose to sketch the support of the probability measure (that does not need to be a manifold) by some gradient flow complex, or more precisely by its Hasse diagram. The gradient flow is defined with respect to the distance function to the sample points. We prove that a gradient flow complex (that can be computed) is homotopy equivalent to the support of the measure for sufficiently dense samplings, and demonstrate the feasibility of our approach on real world data sets.

Algorithmically, many machine learning tasks boil down to solving parameterized optimization problems. Finding good values for the parameters has significant influence on the statistical performance of these methods. Thus supporting the choice of parameter values algorithmically has received quite some attention recently, especially algorithms for computing the whole solution path of parameterized optimization problem. These algorithms can be used, for instance, to track the solution of a regularized learning problem along the regularization parameter path, or for tracking the solution of kernelized problems along a kernel hyperparameter path. Since exact path following algorithms can be numerically unstable, robust and efficient approximate path tracking algorithms became popular for regularized learning problems. By now algorithms with optimal path complexity are known for many regularized learning problems. That is not the case for kernel hyperparameter path tracking algorithms, where the exact path tracking algorithms can also suffer from numerical instabilities. The robust approximation algorithms for regularization path tracking can not be used directly for kernel hyperparameter path tracking problems since the latter fall into a different problem class. Here we address this problem by devising a robust and efficient path tracking algorithm that can also handle kernel hyperparameter paths and has asymptotically optimal complexity. We use this algorithm to compute approximate kernel hyperparamter solution paths for support vector machines and robust kernel regression. Experimental results for this problem applied to various data sets confirms the theoretical complexity analysis.

We present a parallel algorithm and its implementation for computing the entire Hasse diagram of the flow complex of a point cloud in Euclidean space. Known algorithms for computing the flow complex in two and three dimensions compute the geometric realization of the flow complex and need to compute the Delaunay triangulation of the point cloud first. Our algorithm computes less information, namely only the Hasse diagram of the flow complex that is augmented with enough geometric information to allow the same topological multi-scale analysis of point cloud data as the alpha shape filtration without computing the Delaunay triangulation explicitly. We show experimental results for medium dimensions that demonstrate that our algorithm scales well with the number of available cores on a multicore architecture.

We consider an abstract class of optimization problems that are parameterized concavely in a single parameter, and show that the solution path along the parameter can always be approximated with accuracy $\varepsilon >0$ by a set of size $O(1/\sqrt{\varepsilon})$. A lower bound of size $\Omega (1/\sqrt{\varepsilon})$ shows that the upper bound is tight up to a constant factor. We also devise an algorithm that calls a step-size oracle and computes an approximate path of size $O(1/\sqrt{\varepsilon})$. Finally, we provide an implementation of the oracle for soft-margin support vector machines, and a parameterized semi-definite program for matrix completion.

We consider parameterized convex optimization problems over the unit simplex, that depend on one parameter. We provide a simple and efficient scheme for maintaining an ϵ-approximate solution (and a corresponding ϵ-coreset) along the entire parameter path. We prove correctness and optimality of the method. Practically relevant instances of the abstract parameterized optimization problem are for example regularization paths of support vector machines, multiple kernel learning, and minimum enclosing balls of moving points.

Color mapping and semitransparent layering play an important role in many visualization scenarios, such as information visualization and volume rendering. The combination of color and transparency is still dominated by standard alpha-compositing using the Porter-Duff over operator which can result in false colors with deceiving impact on the visualization. Other more advanced methods have also been proposed, but the problem is still far from being solved. Here we present an alternative to these existing methods specifically devised to avoid false colors and preserve visual depth ordering. Our approach is data driven and follows the recently formulated knowledge-assisted visualization (KAV) paradigm. Preference data, that have been gathered in web-based user surveys, are used to train a support-vector machine model for automatically predicting an optimized hue-preserving blending. We have applied the resulting model to both volume rendering and a specific information visualization technique, illustrative parallel coordinate plots. Comparative renderings show a significant improvement over previous approaches in the sense that false colors are completely removed and important properties such as depth ordering and blending vividness are better preserved. Due to the generality of the defined data-driven blending operator, it can be easily integrated also into other visualization frameworks.

We show that conjoint analysis, a popular multi-attribute preference assessment technique used in market research, is a well suited tool to evaluate a multitude of gamut mapping algorithms simultaneously. Our analysis is based on data from psycho-visual tests assessed in a laboratory and in a web environment. Conjoint analysis allows us to quantify the contribution of every single parameter value to the perceived value of the algorithm; it also allows us to test the influence of additional parameters like gamut size or color shifts. We show that conjoint analysis can be individualized to images or observers if enough data is available. Especially promising in this respect is the combination of individual and population data.

This paper addresses the fundamental problem of computing stable medial representations of 3D shapes. We propose a spatially adaptive classification of geometric features that yields a robust algorithm for generating medial representations at different levels of abstraction. The recently introduced continuous scale axis transform serves as the mathematical foundation of our algorithm. We show how geometric and topological properties of the continuous setting carry over to discrete shape representations. Our method combines scaling operations of medial balls for geometric simplification with filtrations of the medial axis and provably good conversion steps to and from union of balls, to enable efficient processing of a wide variety shape representations including polygon meshes, 3D images, implicit surfaces, and point clouds. We demonstrate the robustness and versatility of our algorithm with an extensive validation on hundreds of shapes including complex geometries consisting of millions of triangles.

We demonstrate how the scale axis transform can be used to compute a parameterized family of shape skeletons. The skeletons gradually represent only the most important features of a shape, in a scale-adaptive manner. Here a shape O is any bounded open subset of the plane R2. The scale axis for scale value $s$ is the medial axis of the multiplicatively grown shape O_s, where Os is the union of medial balls of O with radii scaled by the factor s.

We introduce the scale axis transform, a new skeletal shape representation for bounded open sets O ⊂ Rd. The scale axis transform induces a family of skeletons that captures the important features of a shape in a scale-adaptive way and yields a hierarchy of successively simplified skeletons. Its definition is based on the medial axis transform and the simplification of the shape under multiplicative scaling: the s-scaled shape Os is the union of the medial balls of O with radii scaled by a factor of s. The s-scale axis transform of O is the medial axis transform of Os, with radii scaled back by a factor of 1/s. We prove topological properties of the scale axis transform and we describe the evolution s → Os by defining the multiplicative distance function to the shape and studying properties of the corresponding steepest ascent flow. All our theoretical results hold for any dimension. In addition, using a discrete approximation, we present several examples of two-dimensional scale axis transforms that illustrate the practical relevance of our new framework.

Professional designers and artists are quite cognizant of the rules that guide the design of effective color palettes, from both aesthetic and attention-guiding points of view. In the field of visualization, however, the use of systematic rules embracing these aspects has received less attention. The situation is further complicated by the fact that visualization often uses semi-transparencies to reveal occluded objects, in which case the resulting color mixing effects add additional constraints to the choice of the color palette. Color design forms a crucial part in visual aesthetics. Thus, the consideration of these issues can be of great value in the emerging field of illustrative visualization. We describe a knowledge-based system that captures established color design rules into a comprehensive interactive framework, aimed to aid users in the selection of colors for scene objects and incorporating individual preferences, importance functions, and overall scene composition. Our framework also offers new knowledge and solutions for the mixing, ordering and choice of colors in the rendering of semi-transparent layers and surfaces. All design rules are evaluated via user studies, for which we extend the method of conjoint analysis to task-based testing scenarios. Our framework's use of principles rooted in color design with application for the illustration of features in pre-classified data distinguishes it from existing systems which target the exploration of continuous-range density data via perceptual color maps.

We illustrate a simple algorithm for approximating the medial axis of a 2D shape with smooth boundary from a sample of this boundary. The algorithm is compared to a more general approximation method that builds on the same idea, namely, to approximate the shape by a union of balls. While not as general, our algorithm is simpler, faster and numerically more stable. Both algorithms are visualized using the Mesecina tool, which is also described.

We explore the potential of image-dependent gamut mapping as a constrained optimization problem. The performance of our new approach is compared to standard reference gamut mapping algorithms in psycho-visual tests.

Visualization algorithms can have a large number of parameters, making the space of possible rendering results rather high-dimensional. Only a systematic analysis of the perceived quality can truly reveal the optimal setting for each such parameter. However, an exhaustive search in which all possible parameter permutations are presented to each user within a study group would be infeasible to conduct. Additional complications may result from possible parameter co-dependencies. Here, we will introduce an efficient user study design and analysis strategy that is geared to cope with this problem. The user feedback is fast and easy to obtain and does not require exhaustive parameter testing. To enable such a framework we have modified a preference measuring methodology, conjoint analysis, that originated in psychology and is now also widely used in market research. We demonstrate our framework by a study that measures the perceived quality in volume rendering within the context of large parameter spaces.

The medial axis of a shape is known to carry a lot of information about it. In particular a recent result of Lieutier establishes that every bounded open subset of Rn has the same homotopy type as its medial axis. In this paper we provide an algorithm that, given a sufficiently dense but not necessarily uniform sample from the surface of a shape with smooth boundary, computes a core for its medial axis approximation, in form of a piecewise linear cell complex, that captures the topology of the medial axis of the shape. We also provide a natural method to freely augment this core in order to enhance it geometrically all the while maintaining its topological guarantees. The definition of the core and its extension method are based on the steepest ascent flow induced by the distance function to the sample. We also provide a geometric guarantee on the closeness of the core and the actual medial axis.

We show that any comparison based, randomized algorithm to approximate any given ranking of n items within expected Spearman’s footrule distance n 2/ν(n) needs at least n (min{log ν(n), log n} – 6) comparisons in the worst case. This bound is tight up to a constant factor since there exists a deterministic algorithm that shows that 6n(log ν(n)+1) comparisons are always sufficient.

The distance function to surfaces in three dimensions plays a key role in many geometric modeling applications such as medial axis approximations, surface reconstructions, offset computations, feature extractions and others. In most cases, the distance function induced by the surface is approximated by a discrete distance function induced by a discrete sample of the surface. The critical points of the distance function determine the topology of the set inducing the function. However, no earlier theoretical result has linked the critical points of the distance to a sampling of geometric structures to their topological properties. We provide this link by showing that the critical points of the distance function induced by a discrete sample of a surface either lie very close to the surface or near its medial axis and this closeness is quantified with the sampling density. Based on this result, we provide a new flow-complex-based surface reconstruction algorithm that, given a tight ε-sampling of a surface, approximates the surface geometrically, both in Hausdorff distance and normals, and captures its topology.

We describe methods for computing an implicit model of a hypersurface that is given only by a finite sampling. The methods work by mapping the sample points into a reproducing kernel Hilbert space and then deter- mining regions in terms of hyperplanes. 1 Introduction Suppose we are given a finite sampling (in machine learning terms, training data) x1, . . . , xm X , where the domain X is some hypersurface in Euclidean space Rd. The case d = 3 is especially interesting since these days there are many devices, e.g., laser range scanners, that allow the acquisition of point data from the boundary surfaces of solids. For further processing it is often necessary to transform this data into a continu- ous model. Today the most popular approach is to add connectivity information to the data by transforming them into a triangle mesh (see [4] for an example of such a transformation algorithm). But recently also implicit models, where the surface is modeled as the zero set of some sufficiently smooth function, gained some popularity [1]. They bear resemblance to level set methods used in computer vision [6]. One advantage of implicit models is that they easily allow the derivation of higher order differential quantities such as curvatures. Another advantage is that an inside-outside test, i.e., testing whether a query point lies on the bounded or unbounded side of the surface, boils down to determining the sign of a function-evaluation at the query point. Inside-outside tests are important when one wants to intersect two solids. The goal of this paper is, loosely speaking, to find a function which takes the value zero on a surface which (1) contains the training data and (2) is a "reasonable" implicit model of X . To capture properties of its shape even in the above general case, we need to exploit some structure on X . In line with a sizeable amount of recent work on kernel methods [11], we assume that this structure is given by a (positive definite) kernel, i.e., a real valued function Partially supported by the Swiss National Science Foundation under the project "Non-linear manifold learning". Figure 1: In the 2-D toy example depicted, o o the hyperplane w, (x) = separates all o o but one of the points from the origin. The out- . o o lier (x) is associated with a slack variable , o which is penalized in the objective function /||w|| (4). The distance from the outlier to the hy- o w o /||w|| perplane is / w ; the distance between hy- o x ( ) perplane and origin is / w . The latter im- plies that a small w corresponds to a large margin of separation from the origin. k on X X which can be expressed as k(x, x ) = (x), (x ) (1) for some map into a Hilbert space H. The space H is the reproducing kernel Hilbert space (RKHS) associated with k, and is called its feature map. A popular example, in the case where X is a normed space, is the Gaussian (where > 0) x - x 2 k(x, x ) = exp - . (2) 2 2 The advantage of using a positive definite kernel as a similarity measure is that it allows us to construct geometric algorithms in Hilbert spaces. 2 Single-Class SVMs Single-class SVMs were introduced [8, 10] to estimate quantiles C {x X |f (x) [, [} of an unknown distribution P on X using kernel expansions. Here, f (x) = ik(xi, x) - , (3) i where x1, . . . , xm X are unlabeled data generated i.i.d. according to P . The single-class SVM approximately computes the smallest set C C containing a specified fraction of all training examples, where smallness is measured in terms of the norm in the RKHS H associated with k, and C is the family of sets corresponding to half-spaces in H. Depending on the kernel, this notion of smallness will coincide with the intuitive idea that the quantile estimate should not only contain a specified fraction of the training points, but it should also be sufficiently smooth so that the same is approximately true for previously unseen points sampled from P . Let us briefly describe the main ideas of the approach. The training points are mapped into H using the feature map associated with k, and then it is attempted to separate them from the origin with a large margin by solving the following quadratic program: for (0, 1],1 1 1 minimize w 2 + i - (4) wH,R 2 m m ,R i subject to w, (xi) - i, i 0. (5) Since non-zero slack variables i are penalized in the objective function, we can expect that if w and solve this problem, then the decision function, f (x) = sgn ( w, (x) - ) will 1Here and below, bold face greek character denote vectors, e.g., = (1, . . . , m) , and indices i, j by default run over 1, . . . , m. Figure 2: Models computed with a single class SVM using a Gaussian kernel (2). The three examples differ in the value chosen for in the kernel - a large value (0.224 times the diameter of the hemisphere) in the left figure and a small value (0.062 times the diameter of the hemisphere) in the middle and right figure. In the right figure also non-zero slack variables (outliers) were allowed. Note that that the outliers in the right figure correspond to a sharp feature (non-smoothness) in the original surface. equal 1 for most examples xi contained in the training set,2 while the regularization term w will still be small. For an illustration, see Figure 1. The trade-off between these two goals is controlled by a parameter . One can show that the solution takes the form f (x) = sgn ik(xi, x) - , (6) i where the i are computed by solving the dual problem, 1 minimize ij k(xi, xj ) (7) Rm 2 ij 1 subject to 0 i and m i = 1. (8) i Note that according to (8), the training examples contribute with nonnegative weights i 0 to the solution (6). One can show that asymptotically, a fraction of all training examples will have strictly positive weights, and the rest will be zero (the "-property"). In our application we are not primarily interested in a decision function itself but in the boundaries of the regions in input space defined by the decision function. That is, we are interested in f -1(0), where f is the kernel expansion (3) and the points x1, . . . , xm X are sampled from some unknown hypersurface X Rd. We want to consider f -1(0) as a model for X . In the following we focus on the case d = 3. If we assume that the xi are sampled without noise from X which for example is a reasonable assumption for data obtained with a state of the art 3d laser scanning device we should set the slack variables in (4) and (5) to zero. In the dual problem this results in removing the upper constraints on the i in (8). Note that sample points with non-zero slack variable cannot be contained in f -1(0). But also sample points whose image in feature space lies above the optimal hyperplane are not contained in f -1(0) (see Figure 1) -- we will address this in the next section. It turns out that it is useful in practice to allow non-zero slack variables, because they prevent f -1(0) from decomposing into many connected components (see Figure 2 for an illustration). In our experience, one can ensure that the images of all sample points in feature space lie close to (or on) the optimal hyperplane can be achieved by choosing in the Gaussian 2We use the convention that sgn (z) equals 1 for z 0 and -1 otherwise. Figure 3: Two parallel hy- x ( * )o perplanes w, (x) = + o o /* () enclosing all but two . ||w|| o of the points. The outlier o (x()) is associated with (+ * ||w|| )/ (+)/||w|| o a slack variable (), which w o /||w|| o x ( ) o is penalized in the objective function (9). kernel (2) such that the Gaussians in the kernel expansion (3) are highly localized. How- ever, highly localized Gaussians are not well suited for interpolation -- the implicit surface decomposes into several components. Allowing outliers mitigates the situation to a certain extent. Another way to deal with the problem is to further restrict the optimal region in feature space. In the following we will pursue the latter approach.

We introduce the adaptive neighborhood graph as a data structure for modeling a smooth manifold M embedded in some (potentially very high-dimensional) Euclidean space Rd. We assume that M is known to us only through a finite sample P? M, as it is often the case in applications. The adaptive neighborhood graph is a geometric graph on P. Its complexity is at most min[2O(k)n, n2], where n=|P| and k=dim M, as opposed to the n[d/2] complexity of the Delaunay triangulation, which is often used to model manifolds. We show that we can provably correctly infer the connectivity of M and the dimension of M from the adaptive neighborhood graph provided a certain standard sampling condition is fulfilled. The running time of the dimension detection algorithm is d2O(k7log k) for each connected component of M. If the dimension is considered constant, this is a constant-time operation, and the adaptive neighborhood graph is of linear size. Moreover, the exponential dependence of the constants is only on theintrinsic dimension k, not on the ambient dimension d. This is of particular interest if the co-dimension is high, i.e., if k is much smaller than d, as is the case in many applications. The adaptive neighborhood graph also allows us to approximate the geodesic distances between the points in P.

In this paper we establish a topological similarity between two apparently different shape constructors from a set of points. Shape constructors are geometric structures that transform finite point sets into continuous shapes. Due to their immense practical importance in geometric modeling various shape constructors have been proposed recently. Understanding the relations among them often leads to new insights that are potentially helpful in applications. Here we discover a topological equivalence among two such geometric structures, namely α shapes and flow shapes. Both shapes found applications in surface reconstruction and molecular modelin.

A major challenge in requirements engineering is to analyze the stakeholders preferences. Software projects are usually described by a large set of attributes. By varying the values of these attributes, the project could potentially be realized in various different ways. The attributes have to be chosen in order to meet the stakeholders' expectations. Therefore one needs to know the utility of a given project realization as perceived by the stakeholders. Here we introduce a well established technique from marketing research to obtain individual utility functions for each stakeholder. These functions also provide a natural way to deal with attribute interdependencies.

We present an efficient algorithm that computes a manifold triangular mesh from a set of unorganized sample points in . The algorithm builds on the observation made by several researchers that the Gabriel graph of the sample points provides a good surface description. However, this surface description is only one-dimensional. We associate the edges of the Gabriel graph with index 1 critical points of a dynamical system induced by the sample points. Exploiting also the information contained in the critical points of index 2 provides a two-dimensional surface description which can be easily turned into a manifold.

Current surface reconstruction algorithms perform satisfactorily on we ll-sampled, smooth surfaces without boundaries. However, these algorithms face difficulty with undersampling. Cases of undersampling are prevalent in real data since often they sample a part of the boundary of an object, or are derived from a surface with high curvature or nonsmoothness. In this paper we present an algorithm to detect the boundaries where dense sampling stops and undersampling begins. This information can be used to reconstruct surfaces with boundaries, and also to localize small and sharp features where usually undersampling happens. We report the effectiveness of the algorithm with a number of experimental results. Theoretically, we justify the algorithm with some mild assumptions that are valid for most practical data.

Shape modeling is an integral part of many visualization problems. Recent advances in scanning technology and a number of surface reconstruction algorithms have opened up a new paradigm for modeling shapes from samples. Many of the problems currently faced in this modeling paradigm can be traced back to two anomalies in sampling, namely undersampling and oversampling. Boundaries, non-smoothness and small features create undersampling problems, whereas oversampling leads to too many triangles. We use Voronoi cell geometry as a unified guide to detect undersampling and oversampling. We apply these detections in surface reconstruction and model simplification. Guarantees of the algorithms can be proved. The authors show the success of the algorithms empirically on a number of interesting data sets.

We present a novel approach to surface reconstruction based on the Delaunay complex. First we give a simple and fast algorithm that picks locally a surface at each vertex. For that, we introduce the concept of /spl lambda/-intervals. It turns out that for smooth regions of the surface this method works very well and at difficult parts of the surface yields an output well-suited for postprocessing. As a postprocessing step we propose a topological clean up and a new technique based on linear programming in order to establish a topologically correct surface. These techniques should be useful also for many other reconstruction schemes.

We give necessary and sufficient regularity conditions under which the curve reconstruction problem is solved by a Traveling Salesman path.