Search

\indent Over the last two decades, randomization has emerged as a leading tool for pursuing efficiency in numerical linear algebra. Its benefits can be explained in part by smoothed analysis, where an algorithm that fails spectacularly in certain cases may be unlikely to do so on random -- or randomly perturbed -- inputs. This observation implies a simple framework for developing accurate and efficient randomized algorithms:\ apply a random perturbation and run an existing method, whose worst-case error (or run-time) can be avoided with high probability. \\\indent Recent work of Banks, Garza-Vargas, Kulkarni, and Srivastava applied this framework to the standard eigenvalue problem, developing a randomized algorithm that (with high probability) diagonalizes a matrix in nearly matrix multiplication time [Foundations of Computational Math 2022]. Central to their work is the phenomenon of \textit{pseudospectral shattering}, in which a small Gaussian perturbation regularizes the pseudospectrum of a matrix, with high probability breaking it into disjoint components and allowing classical, divide-and-conquer eigensolvers to run successfully. Prior to their work, no way of accessing the benefits of divide-and-conquer’s natural parallelization was known in general. \\ \indent In this thesis, we extend the work of Banks et al.\ to the generalized eigenvalue problem -- e.g., $Av = \lambda Bv$ for matrices $A,B \in {\mathbb C}^{n \times n}$. Our main contributions can be summarized as follows. \begin{enumerate}[itemsep=0em] \item First, we show that pseudospectral shattering generalizes directly:\ randomly perturbing $A$ and $B$ has a similar regularizing effect on the pseudospectra of the corresponding matrix pencil $(A,B)$ with high probability. \item Building on pseudospectral shattering, we construct and analyze a fast, randomized, divide-and-conquer algorithm for diagonalizing $(A,B)$, which begins by randomly perturbing the inputs. \item Finally, we demonstrate that both

Radiation therapy is a clinical field in which computer simulations play a crucial role. Before patients undergo radiation therapy, an individual treatment plan for each patient needs to be created based on the specifics of their case (a process often referred to as treatment planning). The main aspect of this is setting the control parameters for the treatment machine so that the radioactive dose delivered to the patient is as concentrated to the tumor volume as possible. The inverse problem of determining such appropriate control parameters is typically formulated as an optimization problem, which, considering the complexity of modern treatment machines, requires computerized algorithms to solve accurately. Solving this optimization problem can be a key computational bottleneck in the treatment planning workflow. In many cases, finding a suitable treatment plan is a trial-and-error process, requiring multiple solutions of the optimization problems with different weights and parameters. This thesis proposes different methods to enable the use of high performance computing (HPC) hardware to accelerate the optimization process in radiation therapy. We deal with two main computational aspects of the optimization workflow, the calculation of dose, gradients and objective functions; as well as the optimization solver itself. For dose calculation during optimization, we propose a CUDA kernel for sparse matrix-vector products tuned to dose deposition matrices from proton therapy. For the evaluation of the objective function -- which is often constructed as a weighted sum -- we develop a method to distribute the calculation of the objective function and its gradient across computational nodes using message passing. For the optimization solver itself we first propose a task-based parallel implementation for Cholesky factorization of banded matrices. This can be an important kernel in interior point methods (IPM) for optimization, depending on the structure of the optimizati, Strålterapi är en klinisk gren där datorsimuleringar spelar en viktig roll. Innan patienters strålbehandlingar påbörjas behövs individuella behandlingsplaner skapas baserat på varje patientfalls specifika omständigheter (en process som ofta kallas dosplanering (eng. treament planning)). Huvudaspekten av denna process är att bestämma styrparametrar för behandlingsmaskinen så att stråldosen som levereras till patienten är så koncentrerad till tumören som möjligt. Det inversa problemet att bestämma sådana styrparametrar formuleras typiskt som ett optimeringsproblem, vilket, givet hur komplexa moderna behandlingsmaskiner är, kräver datoralgoritmer för att lösas precist. Att lösa detta optimeringsproblem kan vara en viktig flaskhals i dosplaneringsprocessen. I många fall är det en iterativ process att hitta en lämplig behandlingsplan som kräver att man testar och löser optimeringsproblemet med flera olika vikter och parametrar. I denna avhandling studerar vi metoder för att möjliggöra användandet av hårdvara från högprestandaberäkningar för att accelerera optimeringsprocessen i strålterapi. Vi angriper två olika beräkningsaspekter i lösningen av optimeringsproblemet: beräkningen av stråldos, gradienter och målfunktioner, samt lösningen av optimeringsproblemet i sig. För dosberäkningen som krävs under optimeringen utvecklar vi en beräkningsmetod i CUDA för matris-vektor produkter med glesa matriser som är specifikt anpassad för dosmatriser (eng. dose deposition matrices) från protonterapi. För beräkningen av målfunktionen --- som i många fall är formulerad som en viktad summa --- utvecklar vi en metod för att distribuera beräkningen av målfunktionen och dess gradient över flera beräkningsnoder genom message passing. För optimeringslösaren själv så utvecklar vi först en parallel implementation av Cholesky-faktorisering för bandade matriser med parallelisering baserad på tasks. Just den beräkningsalgoritmen kan vara viktig för många inrepunktsmetoder för optimering, baserat, QC 20140423

Radiation therapy is a clinical field in which computer simulations play a crucial role. Before patients undergo radiation therapy, an individual treatment plan for each patient needs to be created based on the specifics of their case (a process often referred to as treatment planning). The main aspect of this is setting the control parameters for the treatment machine so that the radioactive dose delivered to the patient is as concentrated to the tumor volume as possible. The inverse problem of determining such appropriate control parameters is typically formulated as an optimization problem, which, considering the complexity of modern treatment machines, requires computerized algorithms to solve accurately. Solving this optimization problem can be a key computational bottleneck in the treatment planning workflow. In many cases, finding a suitable treatment plan is a trial-and-error process, requiring multiple solutions of the optimization problems with different weights and parameters. This thesis proposes different methods to enable the use of high performance computing (HPC) hardware to accelerate the optimization process in radiation therapy. We deal with two main computational aspects of the optimization workflow, the calculation of dose, gradients and objective functions; as well as the optimization solver itself. For dose calculation during optimization, we propose a CUDA kernel for sparse matrix-vector products tuned to dose deposition matrices from proton therapy. For the evaluation of the objective function -- which is often constructed as a weighted sum -- we develop a method to distribute the calculation of the objective function and its gradient across computational nodes using message passing. For the optimization solver itself we first propose a task-based parallel implementation for Cholesky factorization of banded matrices. This can be an important kernel in interior point methods (IPM) for optimization, depending on the structure of the optimizati, Strålterapi är en klinisk gren där datorsimuleringar spelar en viktig roll. Innan patienters strålbehandlingar påbörjas behövs individuella behandlingsplaner skapas baserat på varje patientfalls specifika omständigheter (en process som ofta kallas dosplanering (eng. treament planning)). Huvudaspekten av denna process är att bestämma styrparametrar för behandlingsmaskinen så att stråldosen som levereras till patienten är så koncentrerad till tumören som möjligt. Det inversa problemet att bestämma sådana styrparametrar formuleras typiskt som ett optimeringsproblem, vilket, givet hur komplexa moderna behandlingsmaskiner är, kräver datoralgoritmer för att lösas precist. Att lösa detta optimeringsproblem kan vara en viktig flaskhals i dosplaneringsprocessen. I många fall är det en iterativ process att hitta en lämplig behandlingsplan som kräver att man testar och löser optimeringsproblemet med flera olika vikter och parametrar. I denna avhandling studerar vi metoder för att möjliggöra användandet av hårdvara från högprestandaberäkningar för att accelerera optimeringsprocessen i strålterapi. Vi angriper två olika beräkningsaspekter i lösningen av optimeringsproblemet: beräkningen av stråldos, gradienter och målfunktioner, samt lösningen av optimeringsproblemet i sig. För dosberäkningen som krävs under optimeringen utvecklar vi en beräkningsmetod i CUDA för matris-vektor produkter med glesa matriser som är specifikt anpassad för dosmatriser (eng. dose deposition matrices) från protonterapi. För beräkningen av målfunktionen --- som i många fall är formulerad som en viktad summa --- utvecklar vi en metod för att distribuera beräkningen av målfunktionen och dess gradient över flera beräkningsnoder genom message passing. För optimeringslösaren själv så utvecklar vi först en parallel implementation av Cholesky-faktorisering för bandade matriser med parallelisering baserad på tasks. Just den beräkningsalgoritmen kan vara viktig för många inrepunktsmetoder för optimering, baserat, QC 20140423

Radiation therapy is a clinical field in which computer simulations play a crucial role. Before patients undergo radiation therapy, an individual treatment plan for each patient needs to be created based on the specifics of their case (a process often referred to as treatment planning). The main aspect of this is setting the control parameters for the treatment machine so that the radioactive dose delivered to the patient is as concentrated to the tumor volume as possible. The inverse problem of determining such appropriate control parameters is typically formulated as an optimization problem, which, considering the complexity of modern treatment machines, requires computerized algorithms to solve accurately. Solving this optimization problem can be a key computational bottleneck in the treatment planning workflow. In many cases, finding a suitable treatment plan is a trial-and-error process, requiring multiple solutions of the optimization problems with different weights and parameters. This thesis proposes different methods to enable the use of high performance computing (HPC) hardware to accelerate the optimization process in radiation therapy. We deal with two main computational aspects of the optimization workflow, the calculation of dose, gradients and objective functions; as well as the optimization solver itself. For dose calculation during optimization, we propose a CUDA kernel for sparse matrix-vector products tuned to dose deposition matrices from proton therapy. For the evaluation of the objective function -- which is often constructed as a weighted sum -- we develop a method to distribute the calculation of the objective function and its gradient across computational nodes using message passing. For the optimization solver itself we first propose a task-based parallel implementation for Cholesky factorization of banded matrices. This can be an important kernel in interior point methods (IPM) for optimization, depending on the structure of the optimizati, Strålterapi är en klinisk gren där datorsimuleringar spelar en viktig roll. Innan patienters strålbehandlingar påbörjas behövs individuella behandlingsplaner skapas baserat på varje patientfalls specifika omständigheter (en process som ofta kallas dosplanering (eng. treament planning)). Huvudaspekten av denna process är att bestämma styrparametrar för behandlingsmaskinen så att stråldosen som levereras till patienten är så koncentrerad till tumören som möjligt. Det inversa problemet att bestämma sådana styrparametrar formuleras typiskt som ett optimeringsproblem, vilket, givet hur komplexa moderna behandlingsmaskiner är, kräver datoralgoritmer för att lösas precist. Att lösa detta optimeringsproblem kan vara en viktig flaskhals i dosplaneringsprocessen. I många fall är det en iterativ process att hitta en lämplig behandlingsplan som kräver att man testar och löser optimeringsproblemet med flera olika vikter och parametrar. I denna avhandling studerar vi metoder för att möjliggöra användandet av hårdvara från högprestandaberäkningar för att accelerera optimeringsprocessen i strålterapi. Vi angriper två olika beräkningsaspekter i lösningen av optimeringsproblemet: beräkningen av stråldos, gradienter och målfunktioner, samt lösningen av optimeringsproblemet i sig. För dosberäkningen som krävs under optimeringen utvecklar vi en beräkningsmetod i CUDA för matris-vektor produkter med glesa matriser som är specifikt anpassad för dosmatriser (eng. dose deposition matrices) från protonterapi. För beräkningen av målfunktionen --- som i många fall är formulerad som en viktad summa --- utvecklar vi en metod för att distribuera beräkningen av målfunktionen och dess gradient över flera beräkningsnoder genom message passing. För optimeringslösaren själv så utvecklar vi först en parallel implementation av Cholesky-faktorisering för bandade matriser med parallelisering baserad på tasks. Just den beräkningsalgoritmen kan vara viktig för många inrepunktsmetoder för optimering, baserat, QC 20140423

Radiation therapy is a clinical field in which computer simulations play a crucial role. Before patients undergo radiation therapy, an individual treatment plan for each patient needs to be created based on the specifics of their case (a process often referred to as treatment planning). The main aspect of this is setting the control parameters for the treatment machine so that the radioactive dose delivered to the patient is as concentrated to the tumor volume as possible. The inverse problem of determining such appropriate control parameters is typically formulated as an optimization problem, which, considering the complexity of modern treatment machines, requires computerized algorithms to solve accurately. Solving this optimization problem can be a key computational bottleneck in the treatment planning workflow. In many cases, finding a suitable treatment plan is a trial-and-error process, requiring multiple solutions of the optimization problems with different weights and parameters. This thesis proposes different methods to enable the use of high performance computing (HPC) hardware to accelerate the optimization process in radiation therapy. We deal with two main computational aspects of the optimization workflow, the calculation of dose, gradients and objective functions; as well as the optimization solver itself. For dose calculation during optimization, we propose a CUDA kernel for sparse matrix-vector products tuned to dose deposition matrices from proton therapy. For the evaluation of the objective function -- which is often constructed as a weighted sum -- we develop a method to distribute the calculation of the objective function and its gradient across computational nodes using message passing. For the optimization solver itself we first propose a task-based parallel implementation for Cholesky factorization of banded matrices. This can be an important kernel in interior point methods (IPM) for optimization, depending on the structure of the optimizati, Strålterapi är en klinisk gren där datorsimuleringar spelar en viktig roll. Innan patienters strålbehandlingar påbörjas behövs individuella behandlingsplaner skapas baserat på varje patientfalls specifika omständigheter (en process som ofta kallas dosplanering (eng. treament planning)). Huvudaspekten av denna process är att bestämma styrparametrar för behandlingsmaskinen så att stråldosen som levereras till patienten är så koncentrerad till tumören som möjligt. Det inversa problemet att bestämma sådana styrparametrar formuleras typiskt som ett optimeringsproblem, vilket, givet hur komplexa moderna behandlingsmaskiner är, kräver datoralgoritmer för att lösas precist. Att lösa detta optimeringsproblem kan vara en viktig flaskhals i dosplaneringsprocessen. I många fall är det en iterativ process att hitta en lämplig behandlingsplan som kräver att man testar och löser optimeringsproblemet med flera olika vikter och parametrar. I denna avhandling studerar vi metoder för att möjliggöra användandet av hårdvara från högprestandaberäkningar för att accelerera optimeringsprocessen i strålterapi. Vi angriper två olika beräkningsaspekter i lösningen av optimeringsproblemet: beräkningen av stråldos, gradienter och målfunktioner, samt lösningen av optimeringsproblemet i sig. För dosberäkningen som krävs under optimeringen utvecklar vi en beräkningsmetod i CUDA för matris-vektor produkter med glesa matriser som är specifikt anpassad för dosmatriser (eng. dose deposition matrices) från protonterapi. För beräkningen av målfunktionen --- som i många fall är formulerad som en viktad summa --- utvecklar vi en metod för att distribuera beräkningen av målfunktionen och dess gradient över flera beräkningsnoder genom message passing. För optimeringslösaren själv så utvecklar vi först en parallel implementation av Cholesky-faktorisering för bandade matriser med parallelisering baserad på tasks. Just den beräkningsalgoritmen kan vara viktig för många inrepunktsmetoder för optimering, baserat, QC 20140423

In this paper, we extend a recently proposed approach for inverse scattering with Neumann boundary conditions [Druskin et al., Inverse Probl. 37, 075003 (2021)] to the 1D Schrödinger equation with impedance (Robin) boundary conditions. This method approaches inverse scattering in two steps: first, to extract a reduced order model (ROM) directly from the data and, subsequently, to extract the scattering potential from the ROM. We also propose a novel data-assimilation (DA) inversion method based on the ROM approach, thereby avoiding the need for a Lanczos-orthogonalization (LO) step. Furthermore, we present a detailed numerical study and A comparison of the accuracy and stability of the DA and LO methods.

The Kaczmarz method is an iterative method for solving linear systems of equations. The Kaczmarz method has been around since it was developed by Kaczmarz 1937. The main idea behind the original Kaczmarz method is to orthogonally project the previous x_k onto the solution space given by a row of the system. The block Kaczmarz on the other hand orthogonally projects the previous x_k onto the solution space given by a sub system of equations. Both the original Kaczmarz method and block Kaczmarz method can only solve consistent systems, however the extended Kaczmarz method is an adaptation that makes it possible to solve inconsistent systems. We will look at both deterministic and randomized, row and block selection processes then compare them on both consistent and inconsistent systems of equations.

One of the most well-established codes for modeling non-linear Magnetohydrodynamics (MHD) for tokamak reactors is JOREK, which solves these equations with a Bézier surface based finite element method. This code produces a highly sparse but also very large linear system. The main solver behind the code uses the Generalized Minimum Residual Method (GMRES) with a physics-based preconditioner. Even with the preconditioner there are issues with memory and computation costs and the solver doesn’t always converge well. This work contains the first thorough study of the mathematical properties of the underlying linear system, enabling us to diagnose and pinpoint the cause of hampered convergence. In particular, analyzing the spectral properties of the matrix and the preconditioned system with numerical linear algebra techniques will open the door to research and investigate more performant solver strategies, such as projection methods., Computer Simulations for Science and Engineering (COSSE)

Whereas aberrating layers are typically viewed as an impediment to medical ultrasound imaging, they can, surprisingly, also be used to our benefit. As long as we can model the effect of an aberrating layer, we can utilize ‘model-based imaging’, the imaging technique explored throughout this thesis, to reconstruct ultrasound images where traditional beamforming methods would fail, employing the ever increasing computational power available to us nowadays. Not only does this allow us to image through layers, but it also leads to interesting applications, such as 3D ultrasound imaging with spatially undersampled data, using an aberrating ‘coding mask’. The formulation of a measurementmodel, a fundamental part ofmodel-based imaging, also gives insight into the imaging problem mathematically, and allows us to investigate methods for estimating the effect of an aberrating layer ‘blindly’, i.e., without explicitly measuring it. In this thesis, we thus investigate (a), imaging through a layer when the layer’s aberration effect is known, and how it can be applied to imaging with spatially undersampled data, and (b), methods and algorithms for estimating the effect of the aberrating layer without knowing it a priori. In the first part of this thesis, we illustrate how using model-based imaging can be utilized for 3D ultrasound imaging using a single ultrasound transducer, and equipping it with a plastic coding mask. The plastic mask acts as an analog coder, that scrambles the transmitted and received waves in a manner that is location dependent. As a result, the temporal shape of an ultrasound echo can be used instead of the traditional method of using phase differences between sensors in a sensor array. Imaging is instead accomplished using model-based imaging. By measuring the pulse-echo response of each pixel, we can form an image by solving a regularized linear least squares problem, which takes into account the measured pixel-specific pulse-echo signals. The pr, Signal Processing Systems

[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.

[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.

[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.

Recent years have seen rapid progress in the field of quantum algorithms for linear algebra problems. These algorithms can solve important problems in quantum chemistry, condensed matter physics, and quantum field theory simulation, with potentially exponential speedup compared to classical algorithms. The progress is in part due to the development of new methods to implement matrix functions on a quantum computer. Examples include the linear combination of unitaries (LCU) method, quantum signal processing (QSP), and quantum singular value transformation (QSVT), which is closely related to QSP. Using these methods, one can construct matrix functions to filter out unwanted eigenstates of a given Hermitian matrix, and thereby obtain a target eigenstate on a quantum computer. This technique we call quantum eigenstate filtering.We focus our attention on two specific problems: ground state preparation and energy estimation, and solving linear systems. For the ground state preparation problem, we want to prepare the ground state, i.e., the lowest eigenstate of the Hamiltonian $H$, while for ground state energy estimation, we want to estimate the ground state energy, i.e., the the lowest eigenvalue of the Hamiltonian $H$. The ground state energy estimation problem is cannot be solved in polynomial time even on a quantum computer in the worst case. However, this problem becomes tractable with additional assumptions, which may be satisfied in many real-world applications. In this work we assume access to a good initial guess of the ground state, such that its overlap with the ground state is lower bounded by a parameter $\gamma$. Under this assumption, we present quantum algorithms based on the quantum eigenstate filtering technique to estimate the ground state energy and to prepare the ground state, with nearly optimal scaling with respect to both $\gamma$ and the precision.The above mentioned algorithms work in the setting where the quantum computers are fully fault-tolera

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.

This course provides an introduction to numerical methods and computational techniques arising in aerospace engineering. Applications are drawn from aerospace structures, aerodynamics, dynamics and control, and aerospace systems. Techniques covered include numerical integration of systems of ordinary differential equations; numerical discretization of partial differential equations; and probabilistic methods for quantifying the impact of variability. Specific emphasis is given to finite volume methods in fluid mechanics, and finite element methods in structural mechanics.Acknowledgement: Prof. David Darmofal taught this course in prior years, and created some of the materials found in this OCW site.

This course provides an introduction to numerical methods and computational techniques arising in aerospace engineering. Applications are drawn from aerospace structures, aerodynamics, dynamics and control, and aerospace systems. Techniques covered include numerical integration of systems of ordinary differential equations; numerical discretization of partial differential equations; and probabilistic methods for quantifying the impact of variability. Specific emphasis is given to finite volume methods in fluid mechanics, and finite element methods in structural mechanics.Acknowledgement: Prof. David Darmofal taught this course in prior years, and created some of the materials found in this OCW site.

This course provides an introduction to numerical methods and computational techniques arising in aerospace engineering. Applications are drawn from aerospace structures, aerodynamics, dynamics and control, and aerospace systems. Techniques covered include numerical integration of systems of ordinary differential equations; numerical discretization of partial differential equations; and probabilistic methods for quantifying the impact of variability. Specific emphasis is given to finite volume methods in fluid mechanics, and finite element methods in structural mechanics.Acknowledgement: Prof. David Darmofal taught this course in prior years, and created some of the materials found in this OCW site.

We improve the current best running time value to invert sparse matrices over finite fields, lowering it to an expected O(n^{2.2131}) time for the current values of fast rectangular matrix multiplication. We achieve the same running time for the computation of the rank and nullspace of a sparse matrix over a finite field. This improvement relies on two key techniques. First, we adopt the decomposition of an arbitrary matrix into block Krylov and Hankel matrices from Eberly et al. (ISSAC 2007). Second, we show how to recover the explicit inverse of a block Hankel matrix using low displacement rank techniques for structured matrices and fast rectangular matrix multiplication algorithms. We generalize our inversion method to block structured matrices with other displacement operators and strengthen the best known upper bounds for explicit inversion of block Toeplitz-like and block Hankel-like matrices, as well as for explicit inversion of block Vandermonde-like matrices with structured blocks. As a further application, we improve the complexity of several algorithms in topological data analysis and in finite group theory.

This paper presents an adaptive randomized algorithm for computing the butterfly factorization of an m × n matrix with m ≈ n provided that both the matrix and its transpose can be rapidly applied to arbitrary vectors. The resulting factorization is composed of O(log n) sparse factors, each containing O(n) nonzero entries. The factorization can be attained using O(n3/2 log n) computation and O(n log n) memory resources. The proposed algorithm can be implemented in parallel and can apply to matrices with strong or weak admissibility conditions arising from surface integral equation solvers as well as multi-frontal-based finite-difference, finite-element, or finite-volume solvers. A distributed-memory parallel implementation of the algorithm demonstrates excellent scaling behavior.

This paper presents an adaptive randomized algorithm for computing the butterfly factorization of an m × n matrix with m ≈ n provided that both the matrix and its transpose can be rapidly applied to arbitrary vectors. The resulting factorization is composed of O(log n) sparse factors, each containing O(n) nonzero entries. The factorization can be attained using O(n3/2 log n) computation and O(n log n) memory resources. The proposed algorithm can be implemented in parallel and can apply to matrices with strong or weak admissibility conditions arising from surface integral equation solvers as well as multi-frontal-based finite-difference, finite-element, or finite-volume solvers. A distributed-memory parallel implementation of the algorithm demonstrates excellent scaling behavior.

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.

A very common type of problem within mathematics and numerical analysis are fixed-point problems, which can arise as sub-problems of optimization methods, differential equations solvers and much more. The most basic iterative approach for fixed-point problems is fixed-point iteration, special cases of which actually date back as far as the Babylonians, where it was used to to find the square roots of positive numbers. An issue with fixed-point iteration is that it can be very slow, as a consequence, acceleration methods have been developed, which are, not surprisingly, methods for speeding up fixed-point iteration. One of these methods are called Anderson acceleration, which has a very strong relationship with GMRES, an algorithm for solving systems of linear equations, which at first glance, seems completely unrelated. The purpose of this thesis is to investigate the theory behind this relationship and to test it numerically., Two types of problems which often arise within mathematics and numerical analysis are fixed-point problems and systems of linear equations, therefore, developing algorithms that solve these problems efficiently, is of particular interest. Fixed-point iteration arguably the simplest way of solving fixed-point problem and special cases of this method actually date back as far as the Babylonians, where it was used to to find the square roots of positive numbers. An issue with fixed-point iteration is that it often too slow at finding a solution, therefore acceleration methods, such as Anderson acceleration have been developed. What is interesting is that Anderson acceleration has a very strong relation to "the generalized minimal residual method" (GMRES), a very successful algorithm of solving a system of linear equations, which at first glance, seems completely unrelated. In this thesis we investigate the relationship between these algorithms and test it numerically.

A very common type of problem within mathematics and numerical analysis are fixed-point problems, which can arise as sub-problems of optimization methods, differential equations solvers and much more. The most basic iterative approach for fixed-point problems is fixed-point iteration, special cases of which actually date back as far as the Babylonians, where it was used to to find the square roots of positive numbers. An issue with fixed-point iteration is that it can be very slow, as a consequence, acceleration methods have been developed, which are, not surprisingly, methods for speeding up fixed-point iteration. One of these methods are called Anderson acceleration, which has a very strong relationship with GMRES, an algorithm for solving systems of linear equations, which at first glance, seems completely unrelated. The purpose of this thesis is to investigate the theory behind this relationship and to test it numerically., Two types of problems which often arise within mathematics and numerical analysis are fixed-point problems and systems of linear equations, therefore, developing algorithms that solve these problems efficiently, is of particular interest. Fixed-point iteration arguably the simplest way of solving fixed-point problem and special cases of this method actually date back as far as the Babylonians, where it was used to to find the square roots of positive numbers. An issue with fixed-point iteration is that it often too slow at finding a solution, therefore acceleration methods, such as Anderson acceleration have been developed. What is interesting is that Anderson acceleration has a very strong relation to "the generalized minimal residual method" (GMRES), a very successful algorithm of solving a system of linear equations, which at first glance, seems completely unrelated. In this thesis we investigate the relationship between these algorithms and test it numerically.

We show that the spline kernel of order $p$ is a so-called semiseparable function with semiseparability rank $p$. A consequence of this is that kernel matrices generated by the spline kernel are rank structured matrices that can be stored and factorized efficiently. We use this insight to derive new recursive algorithms with linear complexity in the number of knots for various kernel matrix computations. We also discuss applications of these algorithms, including smoothing spline regression, Gaussian process regression, and some related hyperparameter estimation problems.

In this dissertation, we study two problems that have the theme of extracting information from lower dimensional samples. A number of recent works have studied algorithms for entrywise $\ell_p$-\emph{low rank approximation}, namely algorithms which given an $n \times d$ matrix $A$ (with $n \geq d$), output a rank-$k$ matrix $B$ minimizing $\|A-B\|_p^p=\sum_{i,j} |A_{i,j} - B_{i,j}|^p$ when $p > 0$; and $\|A-B\|_0 = \sum_{i,j} [A_{i,j} \neq B_{i,j}]$ for $p=0$, where $\|A-B\|_0$ denotes the number of entries $(i,j)$ for which $A_{i,j} \neq B_{i,j}$.For $p = 1$, this is often considered more robust than the SVD, while for $p = 0$ this corresponds to minimizing the number of disagreements, or robust PCA. This problem is known to be NP-hard for $p \in \{0,1\}$, already for $k = 1$, and while there are polynomial time approximation algorithms, their approximation factor is at best $\poly(k)$. It was left open if there was a polynomial-time approximation scheme (PTAS) for $\ell_p$-approximation for any $p \geq 0$. We show the following:\begin{enumerate} \item On the algorithmic side, for $p \in (0,2)$, we use a technique called \emph{sketching} to give the first $n^{\poly(k/\eps)}$ time $(1+\eps)$-approximation algorithm. For $p = 0$, there are various problem formulations, a common one being the binary setting for which $A\in\{0,1\}^{n\times d}$ and $B = U \cdot V$, where $U\in\{0,1\}^{n \times k}$ and $V\in\{0,1\}^{k \times d}$. For this setting, we obtain an algorithm with time $n \cdot d^{\poly(k/\eps)}$. \item On the hardness front, for $p \in (1,2)$, we show under the Small Set Expansion Hypothesis and Exponential Time Hypothesis (ETH), there is no constant factor approximation algorithm running in time $2^{k^{\delta}}$ for a constant $\delta > 0$, showing an exponential dependence on $k$ is necessary. We also show for finite fields of constant size, under the ETH, that any fixed constant factor approximation algorithm requires $2^{k^{\delta}}$ time for a cons

In this dissertation, we study two problems that have the theme of extracting information from lower dimensional samples. A number of recent works have studied algorithms for entrywise $\ell_p$-\emph{low rank approximation}, namely algorithms which given an $n \times d$ matrix $A$ (with $n \geq d$), output a rank-$k$ matrix $B$ minimizing $\|A-B\|_p^p=\sum_{i,j} |A_{i,j} - B_{i,j}|^p$ when $p > 0$; and $\|A-B\|_0 = \sum_{i,j} [A_{i,j} \neq B_{i,j}]$ for $p=0$, where $\|A-B\|_0$ denotes the number of entries $(i,j)$ for which $A_{i,j} \neq B_{i,j}$.For $p = 1$, this is often considered more robust than the SVD, while for $p = 0$ this corresponds to minimizing the number of disagreements, or robust PCA. This problem is known to be NP-hard for $p \in \{0,1\}$, already for $k = 1$, and while there are polynomial time approximation algorithms, their approximation factor is at best $\poly(k)$. It was left open if there was a polynomial-time approximation scheme (PTAS) for $\ell_p$-approximation for any $p \geq 0$. We show the following:\begin{enumerate} \item On the algorithmic side, for $p \in (0,2)$, we use a technique called \emph{sketching} to give the first $n^{\poly(k/\eps)}$ time $(1+\eps)$-approximation algorithm. For $p = 0$, there are various problem formulations, a common one being the binary setting for which $A\in\{0,1\}^{n\times d}$ and $B = U \cdot V$, where $U\in\{0,1\}^{n \times k}$ and $V\in\{0,1\}^{k \times d}$. For this setting, we obtain an algorithm with time $n \cdot d^{\poly(k/\eps)}$. \item On the hardness front, for $p \in (1,2)$, we show under the Small Set Expansion Hypothesis and Exponential Time Hypothesis (ETH), there is no constant factor approximation algorithm running in time $2^{k^{\delta}}$ for a constant $\delta > 0$, showing an exponential dependence on $k$ is necessary. We also show for finite fields of constant size, under the ETH, that any fixed constant factor approximation algorithm requires $2^{k^{\delta}}$ time for a cons

This course offers an advanced introduction to numerical linear algebra. Topics include direct and iterative methods for linear systems, eigenvalue decompositions and QR/SVD factorizations, stability and accuracy of numerical algorithms, the IEEE floating point standard, sparse and structured matrices, preconditioning, linear algebra software. Problem sets require some knowledge of MATLAB®.

With the rapid development of Quantum Computers (QC) and QC Simulators, there will be an increased demand for functioning Quantum Algorithms in the near future. Some of the most ubiquitously useful algorithms are solvers for linear systems of equations. Since the conception of the Quantum Linear Solver Algorithm (QLSA) by Harrow, Hassidim and Lloyd (HHL) in 2009, many improvements have been made, although a generic implementation for arbitrary matrices and vectors is still not available. In this thesis a variant of the HHL QLSA is studied, and the open challenges are investigated. Solutions for two of the challenges, namely the Eigenvalue Inversion subroutine and the Higher-Order Ancilla Rotation subroutine, are discussed. As part of the thesis project, these subroutines have been implemented in the QX Quantum Computer Simulator, and the subroutines are combined to form a complete Quantum Linear Solver (QLS), with the restraint that the implementation for the vector and Hamiltonian of the matrix must be provided by the user. A proof-of-concept QLS by Cao et al. is also implemented in the QX simulator, and using the implementation of the vector and Hamiltonian of Cao et al. the complete solver is tested. In the process of this thesis, a framework for basic Quantum Arithmetic is built providing three variants of Integer Adders, two variants of Integer Subtracters, one Integer Multiplier and one Integer Divider. In addition, gates not natively available in the QX simulator are implemented, and a number of improvements and extensions of algorithms presented in the literature are given, making the described algorithms function on the QX simulator and extending features., Applied Mathematics | Applied Physics

This paper argues that randomized linear sketching is a natural tool for on-the-fly compression of data matrices that arise from large-scale scientific simulations and data collection. The technical contribution consists in a new algorithm for constructing an accurate low-rank approximation of a matrix from streaming data. This method is accompanied by an a priori analysis that allows the user to set algorithm parameters with confidence and an a posteriori error estimator that allows the user to validate the quality of the reconstructed matrix. In comparison to previous techniques, the new method achieves smaller relative approximation errors and is less sensitive to parameter choices. As concrete applications, the paper outlines how the algorithm can be used to compress a Navier--Stokes simulation and a sea surface temperature dataset.

In this dissertation, we study two problems that have the theme of extracting information from lower dimensional samples. A number of recent works have studied algorithms for entrywise $\ell_p$-\emph{low rank approximation}, namely algorithms which given an $n \times d$ matrix $A$ (with $n \geq d$), output a rank-$k$ matrix $B$ minimizing $\|A-B\|_p^p=\sum_{i,j} |A_{i,j} - B_{i,j}|^p$ when $p > 0$; and $\|A-B\|_0 = \sum_{i,j} [A_{i,j} \neq B_{i,j}]$ for $p=0$, where $\|A-B\|_0$ denotes the number of entries $(i,j)$ for which $A_{i,j} \neq B_{i,j}$.For $p = 1$, this is often considered more robust than the SVD, while for $p = 0$ this corresponds to minimizing the number of disagreements, or robust PCA. This problem is known to be NP-hard for $p \in \{0,1\}$, already for $k = 1$, and while there are polynomial time approximation algorithms, their approximation factor is at best $\poly(k)$. It was left open if there was a polynomial-time approximation scheme (PTAS) for $\ell_p$-approximation for any $p \geq 0$. We show the following:\begin{enumerate} \item On the algorithmic side, for $p \in (0,2)$, we use a technique called \emph{sketching} to give the first $n^{\poly(k/\eps)}$ time $(1+\eps)$-approximation algorithm. For $p = 0$, there are various problem formulations, a common one being the binary setting for which $A\in\{0,1\}^{n\times d}$ and $B = U \cdot V$, where $U\in\{0,1\}^{n \times k}$ and $V\in\{0,1\}^{k \times d}$. For this setting, we obtain an algorithm with time $n \cdot d^{\poly(k/\eps)}$. \item On the hardness front, for $p \in (1,2)$, we show under the Small Set Expansion Hypothesis and Exponential Time Hypothesis (ETH), there is no constant factor approximation algorithm running in time $2^{k^{\delta}}$ for a constant $\delta > 0$, showing an exponential dependence on $k$ is necessary. We also show for finite fields of constant size, under the ETH, that any fixed constant factor approximation algorithm requires $2^{k^{\delta}}$ time for a cons

Probing methods for trace estimation of large, sparse matrices has been studied for several decades. In recent years, there has been some work to extend these techniques to instead estimate the diagonal entries of these systems directly. We extend some analysis of trace estimators to their corresponding diagonal estimators, propose a new class of deterministic diagonal estimators which are well-suited to parallel architectures along with heuristic arguments for the design choices in their construction, and conclude with numerical results on diagonal estimation and ordering problems, demonstrating the strengths of our newly-developed methods alongside existing methods.

This course offers an advanced introduction to numerical linear algebra. Topics include direct and iterative methods for linear systems, eigenvalue decompositions and QR/SVD factorizations, stability and accuracy of numerical algorithms, the IEEE floating point standard, sparse and structured matrices, preconditioning, linear algebra software. Problem sets require some knowledge of MATLAB®.

With the rapid development of Quantum Computers (QC) and QC Simulators, there will be an increased demand for functioning Quantum Algorithms in the near future. Some of the most ubiquitously useful algorithms are solvers for linear systems of equations. Since the conception of the Quantum Linear Solver Algorithm (QLSA) by Harrow, Hassidim and Lloyd (HHL) in 2009, many improvements have been made, although a generic implementation for arbitrary matrices and vectors is still not available. In this thesis a variant of the HHL QLSA is studied, and the open challenges are investigated. Solutions for two of the challenges, namely the Eigenvalue Inversion subroutine and the Higher-Order Ancilla Rotation subroutine, are discussed. As part of the thesis project, these subroutines have been implemented in the QX Quantum Computer Simulator, and the subroutines are combined to form a complete Quantum Linear Solver (QLS), with the restraint that the implementation for the vector and Hamiltonian of the matrix must be provided by the user. A proof-of-concept QLS by Cao et al. is also implemented in the QX simulator, and using the implementation of the vector and Hamiltonian of Cao et al. the complete solver is tested. In the process of this thesis, a framework for basic Quantum Arithmetic is built providing three variants of Integer Adders, two variants of Integer Subtracters, one Integer Multiplier and one Integer Divider. In addition, gates not natively available in the QX simulator are implemented, and a number of improvements and extensions of algorithms presented in the literature are given, making the described algorithms function on the QX simulator and extending features., Applied Mathematics | Applied Physics

This course offers an advanced introduction to numerical linear algebra. Topics include direct and iterative methods for linear systems, eigenvalue decompositions and QR/SVD factorizations, stability and accuracy of numerical algorithms, the IEEE floating point standard, sparse and structured matrices, preconditioning, linear algebra software. Problem sets require some knowledge of MATLAB®.

This course offers an advanced introduction to numerical linear algebra. Topics include direct and iterative methods for linear systems, eigenvalue decompositions and QR/SVD factorizations, stability and accuracy of numerical algorithms, the IEEE floating point standard, sparse and structured matrices, preconditioning, linear algebra software. Problem sets require some knowledge of MATLAB®.

[EN] Practical implementation of the QR algorithm: how every linear algebra software library computes eigenvalues and eigenvectors

This BSc thesis focuses on trying to find approximate solutions to partial differential equations using the Fourier collocation method. This method uses Fourier basis functions to approximate the solution to a partial differential equation with periodic boundary conditions. Using Fourier basis functions, one does not need to use large matrices, which makes all computations relatively fast. Another benefit is that for smooth enough initial functions, the error converges very fast. We review some theory of Fourier basis functions, some theory of circulant matrices, and the theory underlying the implementation of the Fourier collocation method. We also describe how one can improve the error of the solution by making some extra calculations before applying the Fourier collocation method, and describe three methods of time stepping. In this BSc thesis, the advection-diffusion equation is used for testing the method, as it can be explicitly solved. We finally present some numerical results. These results confirm in large parts what the theory predicts. However, for some initial functions, the error does not converge as one would expect., Partiella differentialekvationer är ekvationer som beskriver hur tillstånd förändrar sig. Dessa uppstår naturligt när man modellerar omvärlden, från hur molekyler sprider sig i celler till hur galaxer roterar. Att lösa dessa ekvationer är ofta näst intill omöjligt, och man behöver därför approximera lösningar med hjälp av datorer. Det här arbetet fokuserar på en speciell metod för att approximera lösningar, Fourier kollokationsmetoden. Denna metoden använder sig av periodiska funktioner (funktioner som upprepar sig med jämna mellanrum). Om man använder dessa funktioner blir alla beräkningar förhållandevis snabba, vilket är önskvärt. I detta arbete beskriver vi den matematiska teorin som underligger metoden. Vi beskriver sedan metoden och hur man kan förbättra approximationerna. Till sist presenterar vi resultatet från diverse numeriska experiment. Dessa bekräftar till stor del vad den underliggande teorin förutspår, men ett var förvånansvärt.

[EN] Practical implementation of the QR algorithm: how every linear algebra software library computes eigenvalues and eigenvectors

[EN] Practical implementation of the QR algorithm: how every linear algebra software library computes eigenvalues and eigenvectors

This BSc thesis focuses on trying to find approximate solutions to partial differential equations using the Fourier collocation method. This method uses Fourier basis functions to approximate the solution to a partial differential equation with periodic boundary conditions. Using Fourier basis functions, one does not need to use large matrices, which makes all computations relatively fast. Another benefit is that for smooth enough initial functions, the error converges very fast. We review some theory of Fourier basis functions, some theory of circulant matrices, and the theory underlying the implementation of the Fourier collocation method. We also describe how one can improve the error of the solution by making some extra calculations before applying the Fourier collocation method, and describe three methods of time stepping. In this BSc thesis, the advection-diffusion equation is used for testing the method, as it can be explicitly solved. We finally present some numerical results. These results confirm in large parts what the theory predicts. However, for some initial functions, the error does not converge as one would expect., Partiella differentialekvationer är ekvationer som beskriver hur tillstånd förändrar sig. Dessa uppstår naturligt när man modellerar omvärlden, från hur molekyler sprider sig i celler till hur galaxer roterar. Att lösa dessa ekvationer är ofta näst intill omöjligt, och man behöver därför approximera lösningar med hjälp av datorer. Det här arbetet fokuserar på en speciell metod för att approximera lösningar, Fourier kollokationsmetoden. Denna metoden använder sig av periodiska funktioner (funktioner som upprepar sig med jämna mellanrum). Om man använder dessa funktioner blir alla beräkningar förhållandevis snabba, vilket är önskvärt. I detta arbete beskriver vi den matematiska teorin som underligger metoden. Vi beskriver sedan metoden och hur man kan förbättra approximationerna. Till sist presenterar vi resultatet från diverse numeriska experiment. Dessa bekräftar till stor del vad den underliggande teorin förutspår, men ett var förvånansvärt.

Osborne's iteration is a method for balancing n x n matrices which is widely used in linear algebra packages, as balancing preserves eigenvalues and stabilizes their numeral computation. The iteration can be implemented in any norm over R^n, but it is normally used in the L_2 norm. The choice of norm not only affects the desired balance condition, but also defines the iterated balancing step itself. In this paper we focus on Osborne's iteration in any L_p norm, where p < infty. We design a specific implementation of Osborne's iteration in any L_p norm that converges to a strictly epsilon-balanced matrix in O~(epsilon^{-2}n^{9} K) iterations, where K measures, roughly, the number of bits required to represent the entries of the input matrix. This is the first result that proves a variant of Osborne's iteration in the L_2 norm (or any L_p norm, p < infty) strictly balances matrices in polynomial time. This is a substantial improvement over our recent result (in SODA 2017) that showed weak balancing in L_p norms. Previously, Schulman and Sinclair (STOC 2015) showed strict balancing of another variant of Osborne's iteration in the L_infty norm. Their result does not imply any bounds on strict balancing in other norms.

Fast grasping of unknown objects has crucial impact on the efficiency of robot manipulation especially subjected to unfamiliar environments. In order to accelerate grasping speed of unknown objects, principal component analysis is utilized to direct the grasping process. In particular, a single-view partial point cloud is constructed and grasp candidates are allocated along the principal axis. Force balance optimization is employed to analyze possible graspable areas. The obtained graspable area with the minimal resultant force is the best zone for the final grasping execution. It is shown that an unknown object can be more quickly grasped provided that the component analysis principle axis is determined using single-view partial point cloud. To cope with the grasp uncertainty, robot motion is assisted to obtain a new viewpoint. Virtual exploration and experimental tests are carried out to verify this fast gasping algorithm. Both simulation and experimental tests demonstrated excellent performances based on the results of grasping a series of unknown objects. To minimize the grasping uncertainty, the merits of the robot hardware with two 3D cameras can be utilized to suffice the partial point cloud. As a result of utilizing the robot hardware, the grasping reliance is highly enhanced. Therefore, this research demonstrates practical significance for increasing grasping speed and thus increasing robot efficiency under unpredictable environments., Robot Dynamics, Transport Engineering and Logistics

Fast grasping of unknown objects has crucial impact on the efficiency of robot manipulation especially subjected to unfamiliar environments. In order to accelerate grasping speed of unknown objects, principal component analysis is utilized to direct the grasping process. In particular, a single-view partial point cloud is constructed and grasp candidates are allocated along the principal axis. Force balance optimization is employed to analyze possible graspable areas. The obtained graspable area with the minimal resultant force is the best zone for the final grasping execution. It is shown that an unknown object can be more quickly grasped provided that the component analysis principle axis is determined using single-view partial point cloud. To cope with the grasp uncertainty, robot motion is assisted to obtain a new viewpoint. Virtual exploration and experimental tests are carried out to verify this fast gasping algorithm. Both simulation and experimental tests demonstrated excellent performances based on the results of grasping a series of unknown objects. To minimize the grasping uncertainty, the merits of the robot hardware with two 3D cameras can be utilized to suffice the partial point cloud. As a result of utilizing the robot hardware, the grasping reliance is highly enhanced. Therefore, this research demonstrates practical significance for increasing grasping speed and thus increasing robot efficiency under unpredictable environments., Robot Dynamics, Transport Engineering and Logistics

Tensors, linear-algebraic extensions of matrices in arbitrary dimensions, have numerous applications in computer science and computational science. Many tensors are sparse, containing more than 90% zero entries. Efficient algorithms can leverage sparsity to do less work, but the irregular locations of the nonzero entries pose challenges to performance engineers. Many tensor operations such as tensor-vector multiplications can be sped up substantially by breaking the tensor into equally sized blocks (only storing blocks which contain nonzeros) and performing operations in each block using carefully tuned code. However, selecting the best block size is computationally challenging. Previously, Vuduc et al. defined the fill of a sparse tensor to be the number of stored entries in the blocked format divided by the number of nonzero entries, and showed that the fill can be used as an effective heuristic to choose a good block size. However, they gave no accuracy bounds for their method for estimating the fill, and it is vulnerable to adversarial examples. In this paper, we present a sampling-based method for finding a (1 + epsilon)-approximation to the fill of an order N tensor for all block sizes less than B, with probability at least 1 - delta, using O(B^(2N) log(B^N / delta) / epsilon^2) samples for each block size. We introduce an efficient routine to sample for all B^N block sizes at once in O(N B^N) time. We extend our concentration bounds to a more efficient bound based on sampling without replacement, using the recent Hoeffding-Serfling inequality. We then implement our algorithm and compare our scheme to that of Vuduc, as implemented in the Optimized Sparse Kernel Interface (OSKI) library. We find that our algorithm provides faster estimates of the fill at all accuracy levels, providing evidence that this is both a theoretical and practical improvement. Our code is available under the BSD 3-clause license at https://github.com/peterahrens/FillEstimation.

Tensors, linear-algebraic extensions of matrices in arbitrary dimensions, have numerous applications in computer science and computational science. Many tensors are sparse, containing more than 90% zero entries. Efficient algorithms can leverage sparsity to do less work, but the irregular locations of the nonzero entries pose challenges to performance engineers. Many tensor operations such as tensor-vector multiplications can be sped up substantially by breaking the tensor into equally sized blocks (only storing blocks which contain nonzeros) and performing operations in each block using carefully tuned code. However, selecting the best block size is computationally challenging. Previously, Vuduc et al. defined the fill of a sparse tensor to be the number of stored entries in the blocked format divided by the number of nonzero entries, and showed that the fill can be used as an effective heuristic to choose a good block size. However, they gave no accuracy bounds for their method for estimating the fill, and it is vulnerable to adversarial examples. In this paper, we present a sampling-based method for finding a (1 + epsilon)-approximation to the fill of an order N tensor for all block sizes less than B, with probability at least 1 - delta, using O(B^(2N) log(B^N / delta) / epsilon^2) samples for each block size. We introduce an efficient routine to sample for all B^N block sizes at once in O(N B^N) time. We extend our concentration bounds to a more efficient bound based on sampling without replacement, using the recent Hoeffding-Serfling inequality. We then implement our algorithm and compare our scheme to that of Vuduc, as implemented in the Optimized Sparse Kernel Interface (OSKI) library. We find that our algorithm provides faster estimates of the fill at all accuracy levels, providing evidence that this is both a theoretical and practical improvement. Our code is available under the BSD 3-clause license at https://github.com/peterahrens/FillEstimation.

Sparse matrix-matrix multiplication (or SpGEMM) is a key primitive for many high-performance graph algorithms as well as for some linear solvers, such as algebraic multigrid. The scaling of existing parallel implementations of SpGEMM is heavily bound by communication. Even though 3D (or 2.5D) algorithms have been proposed and theoretically analyzed in the at MPI model on Erd}os{Rffenyi matrices, those algorithms had not been implemented in practice and their complexities had not been analyzed for the general case. In this work, we present the first implementation of the 3D SpGEMM formulation that exploits multiple (intranode and internode) levels of parallelism, achieving significant speedups over the state-of-the-art publicly available codes at all levels of concurrencies. We extensively evaluate our implementation and identify bottlenecks that should be subject to further research.

Sparse matrix-matrix multiplication (or SpGEMM) is a key primitive for many high-performance graph algorithms as well as for some linear solvers, such as algebraic multigrid. The scaling of existing parallel implementations of SpGEMM is heavily bound by communication. Even though 3D (or 2.5D) algorithms have been proposed and theoretically analyzed in the at MPI model on Erd}os{Rffenyi matrices, those algorithms had not been implemented in practice and their complexities had not been analyzed for the general case. In this work, we present the first implementation of the 3D SpGEMM formulation that exploits multiple (intranode and internode) levels of parallelism, achieving significant speedups over the state-of-the-art publicly available codes at all levels of concurrencies. We extensively evaluate our implementation and identify bottlenecks that should be subject to further research.

In recent years, quasi-Newton schemes have proven to be a robust and efficient way for the coupling of partitioned multi-physics simulations in particular for fluid-structure interaction. The focus of this work is put on the coupling of partitioned fluid-structure interaction, where minimal interface requirements are assumed for the respective field solvers, thus treated as black box solvers. The coupling is done through communication of boundary values between the solvers. In this thesis a new quasi-Newton variant (IQN-IMVJ) based on a multi-vector update is investigated in combination with serial and parallel coupling systems. Due to implicit incorporation of passed information within the Jacobian update it renders the problem dependent parameter of retained previous time steps unnecessary. Besides, a whole range of coupling schemes are categorized and compared comprehensively with respect to robustness, convergence behaviour and complexity. Those coupling algorithms differ in the structure of the coupling, i.,e., serial or parallel execution of the field solvers and the used quasi-Newton methods. A more in-depth analysis for a choice of coupling schemes is conducted for a set of strongly coupled FSI benchmark problems, using the in-house coupling library preCICE. The superior convergence behaviour and robust nature of the IQN-IMVJ method compared to well known state of the art methods such as the IQN-ILS method, is demonstrated here. It is confirmed that the multi-vector method works optimal without the need of tuning problem dependent parameters in advance. Furthermore, it appears to be especially suitable in conjunction with the parallel coupling system, in that it yields fairly similar results for parallel and serial coupling. Although we focus on FSI simulation, the considered coupling schemes are supposed to be equally applicable to various kinds of different volume- or surface-coupled problems., In den letzten Jahren haben sich quasi-Newton Verfahren als robuste und effiziente Methode für die Kopplung partitionierter Multiphysik-Simulationen herausgestellt, insbesondere im Bereich der Kopplung von Fluid-Struktur Interaktion. Der Fokus dieser Arbeit liegt auf der Kopplung partitionierten Fluid-Struktur Interaktion mit minimalen Schnittstellenanforderungen an die betreffenden Fluid- oder Struktur-Löser, die daher als black-box Löser angesehen werden können. Die Kopplung an sich wird durch Kommunikation von Randwerten zwischen den Lösern realisiert. Im Rahmen dieser Arbeit wird eine neue quasi-Newton Methode in Kombination mit einem seriellen und parallelen Kopplungssystem untersucht. Diese Variante ist im Gegensatz zu bisherigen Methoden weitestgehend frei von problemabhängigen Parametern wie z. B. der optimalen Menge an wiederverwendeter Information aus vergangenen Zeitschritten, indem alle bisher bekannte Information in einer norm-minimalen und impliziten Art und Weise für die Aufdatierung der Jacobi Matrix verwendet wird. Darüberhinaus betrachten wir in dieser Arbeit ein ganzes Spektrum an Kopplungs-Schemata, die verschiedene quasi-Newton Varianten mit seriellen und parallellen Kopplungsansätzen kombinieren. Eine sorgfältige Klassifizierung sowie ein umfassender Vergleich der Verfahren bezüglich Robustheit, Komplexität und des Konvergenzverhaltens, verschaffen einen guten Überblick. Für eine Auswahl der besten Schemata wird eine eigehendere Analyse anhand einer Reihe anspruchsvoller FSI Anwendungen mit Hilfe der Kopplungs-Bibliothek preCICE durchgeführt. In diesem Zusammenhang gehen wir auf die robuste sowie effiziente Implementierung der Kopplungs-Algorithmen ein. Im Zuge der Experimente zeigt die IQN-IMVJ quasi-Newton Methode ein überlegenes und weitaus robusteres Konvergenzverhalten im Vergleich zu bisherigen Varianten wie beispielsweise die IQN-ILS Methode und arbeitet optimal ohne vorherige Justierung problemabhängiger Parameter. Darüber hinaus ist si