Search

We show that the H2-norm model reduction problem for single-input/single-output (SISO) linear time-invariant (LTI) systems is essentially an eigenvalue problem (EP), from which the globally optimal solution(s) can be retrieved. The first-order optimality conditions of this model reduction problem constitute a system of multivariate polynomial equations that can be converted to an affine (or inhomogeneous) multiparameter eigenvalue problem (AMEP). We solve this AMEP by using the so-called augmented block Macaulay matrix, which is introduced in this paper and has a special (block) multi-shift invariant null space. The set of all stationary points of the optimization problem, i.e., the (2r)-tuples (r is the order of the reduced model) of affine eigenvalues and eigenvectors of the AMEP, follows from a standard EP related to the structure of that null space. At least one of these (2r)-tuples corresponds to the globally optimal solution of the H2-norm model reduction problem. We present a simple numerical example to illustrate our approach.

In this work we develop a non-parametric method to estimate the input-passivity index of an unknown linear and time-invariant (LTI) system from iterative experiments based on the power method from numerical linear algebra. Inspired by the power method for estimating the H∞-norm (or L2-gain) from data, we propose an algorithm that time-reverses input-output data in order to emulate measurements of a virtual system whose L2-gain matches the passivity index of the original system under study. While the proposed method requires exciting the original system twice, we also introduce an improved sampling scheme where only one experiment per iteration is needed.

A kind of algebra, called numerical algebra, is proposed and investigated. As its opponent, non-numerical algebra is also defined. The numeralization and dis-numeralization, which convert non-numerical algebra to numerical algebra and vise versa, are considered. Product structure matrix (PSM) of a finite dimensional algebra is constructed. Using PSM, some fundamental properties of finite dimensional algebras are obtained. Then a necessary and sufficient condition for a numerical algebra to be a field is presented. Finally, the invertibility of Segre (commutative) quaternion and some related properties of matrices over Segre quaternion are investigated.

This paper presents a simple novel IPM matrix computation method, which does not require any additional solutions to calculate homography matrices. This method only requires determining the distance from camera to projecting surface. It’s simplicity allows execution even on very small energy-efficient embedded device processors without use of complex libraries such as OpenCV, that cannot be deployed there. Furthermore, it has strictly analytical solution, which allows avoidance of false optimums that are probable in heuristic methods for optimization. We provide formulas, that are easy transferable to any device and can be calculated using standard math libraries.

Eigenvector computation such as Singular Value Decomposition (SVD) is one of the most fundamental problems in machine learning, optimization and numerical linear algebra. In recent years, many stochastic variance reduction algorithms and randomized coordinate descent algorithms have been developed to efficiently solve the leading eigenvalue problem. By taking full advantage of both variance reduction and randomized coordinate descent techniques, this paper proposes a novel Semi-stochastic Block Coordinate Descent algorithm (SBCD-SVD), which is more suitable than existing algorithms for large-scale leading eigenvalue problems of SVD, and can obtain linear convergence. Unlike existing stochastic variance reduction and randomized coordinate descent methods, our algorithm inherits their advantages. Moreover, we propose a new Asynchronous parallel Semi-stochastic Block Coordinate Descent algorithm (ASBCD-SVD) and one new Asynchronous parallel Sparse approximated Variance Reduction algorithm (ASVR-SVD) for large-scale dense and sparse datasets, respectively. Finally, we prove that both dense and sparse asynchronous parallel variants can converge linearly. Extensive experimental results show that our algorithms attain high parallel speedup and achieve almost the same performance with significantly shorter time, and thus they can be widely used in various practice applications.

Electrical capacitance tomography is a widely used sensor modality for flow imaging in many industrial settings. Adaptive electrical capacitance volume tomography (AECVT) extends the capabilities of traditional ECT by enabling direct volumetric imaging and an improved resolution. Construction of the sensitivity matrix is a necessary step to obtain flow images. This step requires the computation of the electric field inside the sensing domain, which is done via a typical field solver, such as the finite-element method. In this work, we present an efficient and flexible method to construct the sensitivity matrix for AECVT based on individual electrode segment excitations and their judicious combination to form desired matrix elements. We illustrate how the proposed method yields the same sensitivity matrix as the traditional method but at a much lower computational cost. Once all segment contributions are obtained, we also indicate how the proposed method, unlike the traditional approach, can generate the sensitivity matrix on demand for an arbitrary combination of synthetic electrodes and obviating the need for any additional field computations. Finally, we present image reconstruction results for two different experimental scenarios where the mutual capacitance data and the corresponding sensitivity vectors are obtained through the proposed measurement combination scheme.

Sparse hyperspectral unmixing has been a hot topic in recent years. Joint sparsity assumes that each pixel in a small neighborhood of hyperspectral images (HSIs) is composed of the same endmembers, which results in a few nonzero rows in the abundance matrix. Recall that a plethora of unmixing algorithms transform a 3-D HSI into a 2-D matrix with vertical priority. The transformation makes matrix computation easier. It is, however, hard to maintain the horizontal spatial information in HSIs in many cases. To make further use of the spatial information of HSIs, in this article, we propose a bilateral joint-sparse structure for hyperspectral unmixing in an attempt to exploit the local joint sparsity of the abundance matrix in both the vertical and horizontal directions. In particular, we introduce a permutation matrix to realize the bilateral joint-sparse representation and there is no need to construct the matrix explicitly. Moreover, we propose to simultaneously impose the bilateral joint-sparse structure and low rankness on the abundance and develop a new algorithm named bilateral joint-sparse and low-rank unmixing. The proposed algorithm is based on the alternating direction method of multipliers framework and employs a reweighting strategy. The convergence analysis of the proposed algorithm is investigated. Simulated and real-data experiments show the effectiveness of the proposed algorithm.

This article mainly investigates the classic resonant cavity problem with anisotropic and nonconductive media, which is a linear vector Maxwell’s eigenvalue problem. The finite-element method based on edge element of the lowest order and standard linear element is used to solve this type of 3-D closed cavity problem. In order to eliminate spurious zero modes in the numerical simulation, the divergence-free condition supported by Gauss’ law is enforced in a weak sense. After the finite-element discretization, the generalized eigenvalue problem with a linear constraint condition needs to be solved. The penalty method, augmented method, and projection method are applied to solve this difficult problem in numerical linear algebra. The advantages and disadvantages of these three computational methods are also given in this article. Furthermore, we prove that the augmented method is free of spurious modes as long as the anisotropic material is not magnetic lossy. The projection method based on singular value decomposition technique can be used to solve the resonant cavity problem. Moreover, the projection method cannot introduce any spurious modes. At last, several numerical experiments are carried out to verify our theoretical results.

Data preconditioning technique, which reduces the condition number of the problem by a linear transformation of the data matrix, is typically used to accelerate the convergence of the first-order optimization methods for regularized loss minimization. One obvious limitation of the technique is exceedingly expensive of computational cost for the large-scale problems, especially an ocean of samples. In this paper, we have a gradient preconditioning trick and combine it with mini-batch SGD. The proposed gradient preconditioned mini-batch SGD algorithm boosts indeed the convergence with lower computational cost than that of the data preconditioning technique for ridge regression. Concretely, we use recent random projection and linear sketching methods to randomly low rank approximate the data matrix, then we can achieve a appropriate preconditioner through numerical linear algebra. Finally, we apply obtained preconditioner to the gradient to reduce computational cost. The experimental results on both synthetic data and real data sets validate the feasibility and effectiveness of our trick and algorithm.

As widespread development of biometrics, concerns about security and privacy are rapidly increasing. Secure matrix computation is one of the most fundamental and useful operations for statistical a...

In this work, we propose a new set of an intelligent classifier of biomedical color images using Quaternion Discrete Radial Tchebichef Moments (QRTM). This moment has shown very high robustness in terms of color image representation. Thus, we propose a new approach for the fast and accurate reconstruction of multi-level and color images. This approach based on the reconstruction of color images by applying matrix computation. In the second step, we propose a new method for the extraction of Quaternion Discrete Radial Tchebichef Moments invariant (QRTMI). This method based on the representation of the extraction of these invariants. The performance of the invariance of these moments under the three types of geometrical transformations (Rotation, translation, scale) are very important. Finally, we present a new model based on a Multilayer perceptron (MLP) for the classification of biomedical color images. The experimental results show that the QRTMI very powerful compared with radial Legendre - Fourier and Quaternion discrete radial Krawtchouk moments invariant for the biomedical colors images using the BreaKHis Database (Breast Tumor).

Many regions of human movement capturing are commonly used. Still, it includes a complicated capturing method, and the obtained information contains missing information invariably due to the human's body or clothing structure. Recovery of motion that aims to recover from degraded observation and the underlying complete sequence of motion is still a difficult task, because the nonlinear structure and the filming property is integrated into the movements. Machine learning model based two-dimensional matrix computation (MM-TDMC) approach demonstrates promising performance in short-term motion recovery problems. However, the theoretical guarantee for the recovery of nonlinear movement information lacks in the two-dimensional matrix computation model developed for linear information. To overcome this drawback, this study proposes MM-TDMC for human motion and dance recovery. The advantages of the machine learning-based Two-dimensional matrix computation model for human motion and dance recovery shows extensive experimental results and comparisons with auto-conditioned recurrent neural network, multimodal corpus, low-rank matrix completion, and kinect sensors methods.

This survey describes probabilistic algorithms for linear algebraic computations, such as factorizing matrices and solving linear systems. It focuses on techniques that have a proven track record for real-world problems. The paper treats both the theoretical foundations of the subject and practical computational issues.Topics include norm estimation, matrix approximation by sampling, structured and unstructured random embeddings, linear regression problems, low-rank approximation, subspace iteration and Krylov methods, error estimation and adaptivity, interpolatory and CUR factorizations, Nyström approximation of positive semidefinite matrices, single-view (‘streaming’) algorithms, full rank-revealing factorizations, solvers for linear systems, and approximation of kernel matrices that arise in machine learning and in scientific computing.

Graph analytics systems rely on parallel models to realize parallel processing of large-scale graphs, in which the most extensively used models are the bulk synchronous parallel (BSP) and the asynchronous parallel (AP) model. The BSP model is simple and deterministic but has difficulties in expressing some graph processing algorithms, while the AP model is flexible to handle various kinds of algorithms but behaves inefficiently on memory-intensive algorithms due to higher scheduling overhead and internal global barriers. Additionally, although the matrix computation has demonstrated its high-performance characteristic in graph processing, most matrix-based systems only support the BSP model. To improve efficiency of the AP model for graph analytics, this paper proposes a barrierless AP model which uses matrix operations to process graphs and removes the global barrier in each superstep of the AP model to accelerate the execution progress of algorithms. The proposed model is evaluated by comparing it with both AP and BSP-based systems including GraphLab, GRACE, and GraphMat. Experimental results show that our model behaves better than the asynchronous systems, and even outperforms the BSP-based system on some memory-intensive graph algorithms as well as computation-intensive graph algorithms.

The current big data era routinely requires the processing of large-scale data on massive distributed computing clusters. In these applications, data sets are often so large that they cannot be housed in the memory and/or the disk of any one computer. Thus, the data and the processing are typically distributed across multiple nodes. Distributed computation is thus a necessity rather than a luxury. The widespread use of such clusters presents several opportunities and advantages over traditional computing paradigms. However, it also presents newer challenges where coding-theoretic ideas have recently had a significant impact. Large-scale clusters (which can be heterogeneous in nature) suffer from the problem of stragglers, which are slow or failed worker nodes in the system. Thus, the overall speed of a computation is typically dominated by the slowest node in the absence of a sophisticated assignment of tasks to the worker nodes.

This chapter provides a mathematical background to define a mathematical model of cellular systems, and design and optimize wireless communications. It describes the approximation problems and establishes a solution in vector spaces. The chapter reviews some important matrix factorization. Most matrix computation from original data is not easy to calculate in an explicit way. Matrix factorization (or decomposition) allows us to rephrase some matrix forms in such a way that they can be solved more easily. Various matrix decompositions covered include lower upper decomposition, Cholesky decomposition, QR decomposition, and singular value decomposition. The chapter explains the least square algorithm in the context of a regression problem. Regression analysis is a useful tool for forecasting, error reduction or data analysis.

Matrix computations are indispensable tools in science and engineering, while the electronic matrix computations suffered from limited speed. Alternatively, the optical methods offered a high-speed solution. The optical matrix computation is impractical unless it is capable of extending to a large scale, reconfiguring for a general purpose and operating with high efficiency. Here, we report and experimentally demonstrate an optical matrix computing processor based on an integrated linear optical network. The proposed photonic processor is capable of performing fundamental matrix computations including XB = C , AB = X and AX = C , where A , B , C are known matrices, and X is the matrix to be solved. An optical PageRank algorithm is further demonstrated based on the matrix computing processor for the first time. Our demonstration offers an optical method to achieve matrix computations and PageRank algorithm. Meanwhile, it suggests great potential for chip-scale fully programmable matrix computations with self-configuring methods.

Nonlinear conjugate gradient methods are among the most preferable and effortless methods to solve smooth optimization problems. Due to their clarity and low memory requirements, they are more desirable for solving large-scale smooth problems. Conjugate gradient methods make use of gradient and the previous direction information to determine the next search direction, and they require no numerical linear algebra. However, the utility of nonlinear conjugate gradient methods has not been widely employed in solving nonsmooth optimization problems. In this paper, a modified nonlinear conjugate gradient method, which achieves the global convergence property and numerical efficiency, is proposed to solve large-scale nonsmooth convex problems. The new method owns the search direction, which generates sufficient descent property and belongs to a trust region. Under some suitable conditions, the global convergence of the proposed algorithm is analyzed for nonsmooth convex problems. The numerical efficiency of the proposed algorithm is tested and compared with some existing methods on some large-scale nonsmooth academic test problems. The numerical results show that the new algorithm has a very good performance in solving large-scale nonsmooth problems.

Matrix balancing is a preprocessing step in linear algebra computations such as the computation of eigenvalues of a matrix. Such computations are known to be numerically unstable if the matrix is unbalanced, that is the L2 norm of some rows and their corresponding columns are different by orders of magnitude. Given an unbalanced matrix A, the goal of matrix balancing is to find an invertible diagonal matrix D such that DAD^{−1} is balanced or approximately balanced in the sense that every row and its corresponding column have the same norm. In thesis, we study a classic iterative algorithm for matrix balancing due to Osborne (1960). The original algorithm was proposed for balancing rows and columns in the L2 norm, and it works by iterating over balancing a row-column pair in fixed round-robin order. Variants of the algorithm for other norms have been heavily studied and are implemented as standard preconditioners in many numerical linear algebra packages. Despite the popularity of Osborne’s algorithm in practice and extensive research on it there were no polynomial-time upper bound on the running time of this algorithm to explain the excellent performance of this algorithm in practice. Recently (in 2015), Schulman and Sinclair, in a first result of its kind for any norm, analyzed the rate of convergence of a variant of Osborne’s algorithm that uses the L∞ norm and a different order of choosing row-column pairs. In this thesis we study matrix balancing in the L1 norm and other Lp norms. We consider two notions of approximately balancing matrices and refer to them as ε-balancing and strict ε-balancing. As the names suggest strict ε-balancing implies ε-balancing. These notions will be defined in the body of the thesis.We prove upper bounds on the convergence rate of Osborne’s algorithm and some of its vari- ants. We prove fast convergence of different variants of the algorithm to an ε-balanced matrix, and propose a variant that converges to a strictly ε-balanced matrix in polynomial time. These results resolve a problem that has been open since Osborne proposed his algorithm in 1960. The following is a summary of our results for any real matrix A = (a_{ij})_{i,j}:1. We propose a simple greedy variant of Osborne’s algorithm and show that it converges to an ε-balanced matrix in K = O(min{ε^{−2} log w, ε^{−1}n^{3/2} log(w/ε)}) iterations that cost a total of O(m + Kn log n) arithmetic operations over O(n log(w/ε))-bit numbers. Here m is the number of non-zero entries of A, and w = \sum_{i,j} |a_{ij} |/a_{min} with a_{min}= min{|a_{ij} | : a_{ij} ̸= 0}.2. We show that the original round-robin implementation of Osborne’s algorithm converges to an ε-balanced matrix in O(ε^{−2}n^2 log w) iterations totaling O(ε^{−2}mn log w) arithmetic operations over O(n log(w/ε))-bit numbers.3. We devise a random implementation of the iteration and show that it converges to an ε-balanced matrix in O(ε^{−2} log w) iterations using O(m + ε^{−2}n log w) arithmetic operations over O(log(wn/ε))-bit numbers.4. We propose a variant of Osborne’s algorithm and prove that it converges to a strictly ε-balanced matrix in O (ε^{−2}n^9 log(wn/ε) log w/ log n) iterations that requireO (ε^{−2}n^10 log(wn/ε) log w/ log n) arithmetic operations over O(n log w/ε)-bit numbers.5. We demonstrate a lower bound of Ω(1/\sqrt{ε}) on the convergence rate of any implementation of the iteration. Thus, the dependence of our upper bounds on 1/ε is in the right ballpark.All our results are proved for balancing in L1 norm, but we observe, through a known trivial reduction, that our results for L1 balancing apply to any Lp norm for all finite p, at the cost of increasing the number of iterations by only a factor of p (or p2 in some cases).

This thesis concerns the optimization and application of low-rank methods, with a special focus on tensor trains (TTs). In particular, we develop methods for computing TT approximations of a given tensor in a variety of low-rank formats and we show how to solve the tensor completion problem for TTs using Riemannian methods. This is then applied to train a machine learning (ML) estimator based on discretized functions. We also study randomized methods for obtaining low-rank approximations of matrices and tensors. Finally, we consider how such randomized methods can be used to solve general linear matrix and tensor equations.

Various applications in numerical linear algebra and computer science are related to selecting the $$r\times r$$ r × r submatrix of maximum volume contained in a given matrix $$A\in \mathbb R^{n\times n}$$ A ∈ R n × n . We propose a new greedy algorithm of cost $$\mathcal O(n)$$ O ( n ) , for the case A symmetric positive semidefinite (SPSD) and we discuss its extension to related optimization problems such as the maximum ratio of volumes. In the second part of the paper we prove that any SPSD matrix admits a cross approximation built on a principal submatrix whose approximation error is bounded by $$(r+1)$$ ( r + 1 ) times the error of the best rank r approximation in the nuclear norm. In the spirit of recent work by Cortinovis and Kressner we derive some deterministic algorithms, which are capable to retrieve a quasi optimal cross approximation with cost $$\mathcal O(n^3)$$ O ( n 3 ) .

[ES] El auge en la aplicación de redes neuronales profundas (RNPs) en una gran variedad de campos científicos ha propiciado su uso no solo en servidores de cómputo sino también en dispositivos de bajo consumo. Los cálculos que se realizan tanto en entrenamiento como en inferencia de las RNPs se descomponen en núcleos de álgebra lineal y se extraen de bibliotecas especializadas como Intel MKL, BLIS, etc. Sin embargo, la memoria requerida por estas bibliotecas puede exceder de la capacidad máxima de estos pequeños dispositivos. Además, la gran variedad de hardware de bajo consumo hace prácticamente imposible contar con núcleos de cómputo optimizados para cada modelo. Una opción para reducir el coste de generación y mantenimiento de estas bibliotecas es la utilización de generadores de código automáticos como Apache TVM. Estas herramientas permiten desarrollar un solo código común para todos los dispositivos y posteriormente generar el código ensamblador para cada uno. Además, con TVM solamente se debe generar los núcleos necesarios para un modelo de RNP concreto evitando utilizar memoria del dispositivo con funcionalidad que no se utiliza para cada caso. En este proyecto se pretende generar núcleos computacionales para distintas arquitecturas de forma automática utilizando Apache TVM con el objetivo de reproducir las necesidades de los RNPs., [EN] The boom in the application of Deep Neural Networis (DNNs) in a wide variety of scientific fields has led to their use not only in compute-intensive servers but also in low-power devices. Many of the computations performed m RNPs, both in training and inference, are decomposed into linear algebra kernels and extracted from specialised libraries such as Intel Mkl or BLIS. However, the amount of memory required by these libraries can exceed the maximum capacity of these small devices. In addition, the wide variety of low-power hardware makes it virtually impossible to have optimised compute cores for each modeL One option to reduce the cost of generating and maintaming these libraries is the use of automatic code generators such as Apache TVM. These tools allow you to generate a single common code for all devices and then ,generate the assembly code for each device. With TVM, only the cores needed for a specific RNP model must be generated, avoiding the use of device memory with functionality that is not going to be used. This project aims to generate optimised computational kernels for different'architectures automatically using Apache TVM with the objective of reproducing the needs of RNP.

Time-resolved serial femtosecond crystallography (TR-SFX) provides access to protein dynamics on sub-picosecond timescales, and with atomic resolution. Due to the nature of the experiment, these datasets are often highly incomplete and the measured diffracted intensities are affected by partiality. To tackle these issues, one established procedure is that of splitting the data into time bins, and averaging the multiple measurements of equivalent reflections within each bin. This binning and averaging often involve a loss of information. Here, we propose an alternative approach, which we call low-pass spectral analysis (LPSA). In this method, the data are projected onto the subspace defined by a set of trigonometric functions, with frequencies up to a certain cutoff. This approach attenuates undesirable high-frequency features and facilitates retrieving the underlying dynamics. A time-lagged embedding step can be included prior to subspace projection to improve the stability of the results with respect to the parameters involved. Subsequent modal decomposition allows to produce a low-rank description of the system's evolution. Using a synthetic time-evolving model with incomplete and partial observations, we analyze the LPSA results in terms of quality of the retrieved signal, as a function of the parameters involved. We compare the performance of LPSA to that of a range of other sophisticated data analysis techniques. We show that LPSA allows to achieve excellent dynamics reconstruction at modest computational cost. Finally, we demonstrate the superiority of dynamics retrieval by LPSA compared to time binning and merging, which is, to date, the most commonly used method to extract dynamical information from TR-SFX data., Structural Dynamics, 9 (4), ISSN:2329-7778

These lecture notes accompany a refresher course in applied mathematics with a focus on numerical concepts (Part I), numerical linear algebra (Part II), numerical analysis, Fourier series and Fourier transforms (Part III), and differential equations (Part IV). Several numerical projects for group work are provided in Part V. In these projects, the tasks are threefold: mathematical modeling, algorithmic design, and implementation. Therein, it is important to draw interpretations of the obtained results and provide measures (Parts I-IV) how to build confidence into numerical findings such intuition, error analysis, convergence analysis, and comparison to manufactured solutions. Both authors have been jointly teaching over several years this class and bring in a unique mixture of their respective teaching and research fields.

We consider the distributed training of large scale neural networks that serve as PDE (partial differential equation) solvers producing full field outputs. We specifically consider neural solvers for the generalized 3D Poisson equation over megavoxel domains. A scalable framework is presented that integrates two distinct advances. First, we accelerate training a large model via a method analogous to the multigrid technique used in numerical linear algebra. Here, the network is trained using a hierarchy of increasing resolution inputs in sequence, analogous to the `V', `W', `F' and `Half-V' cycles used in multigrid approaches. In conjunction with the multi-grid approach, we implement a distributed deep learning framework which significantly reduces the time to solve. We show scalability of this approach on both GPU (Azure VMs on Cloud) and CPU clusters (PSC Bridges2). This approach is deployed to train a generalized 3D Poisson solver that scales well to predict output full field solutions up to the resolution of 512 X 512 X 512 for a high dimensional family of inputs. This strategy opens up the possibility of fast and scalable training of neural PDE solvers on heterogeneous clusters.