Determination and use of spectral embeddings of large-scale systems by substructuring

A device solves for eigenvalues of a matrix system. The device performs a domain decomposition of a matrix system into non-overlapping subdomains and a reordering of matrices of the matrix system. An interface variable projection subspace associated with interface variables of an adjacency graph of the matrix system is created. The interface variables are related to nodes of the adjacency graph which are connected with nodes located in neighboring partitions. An internal variable projection subspace is created that is associated with internal variables of the adjacency graph of the matrix system, wherein the internal variables are related to nodes of the adjacency graph which are connected only to nodes located in the same partition. A projection matrix is built based on the interface variable projection subspace and the internal variable projection subspace. The device determines eigenvalues that solve a Raleigh-Ritz eigenvalue problem utilizing the projection matrix.

BACKGROUND

Disclosed herein is a system and related method for determining and utilizing eigenvectors (embeddings) of large sparse symmetric matrix pencils, which may occur in a variety of technical applications, including, e.g., spectral clustering and frequency response analysis.

This disclosure considers the determination and application of eigenpairs, consisting of eigenvalues and eigenvectors, using an algebraic domain decomposition scheme to compute a selected number of eigenpairs of large and sparse symmetric matrix pencils. These eigenpairs are computed by means of a Rayleigh-Ritz projection. The projection subspace is divided into two separate subspaces, each one associated with one of the two types of variables: interior, and interface. The projection subspace associated with interface variables is built by computing a few of the eigenvectors and associated first and second derivatives of a zeroth-order approximation of the non-linear matrix-valued interface operator.

The projection subspace associated with interior variables is built independently in each subdomain by exploiting local eigenmodes and matrix resolvent approximations. The sought eigenpairs are then approximated by means of a Rayleigh-Ritz projection. The proposed technique may result in computational resource improvements over other processes, such as the shift-and-invert Lanczos, and Automated MultiLevel Substructuring combined with p-way vertex-based partitionings.

SUMMARY

According to one aspect disclosed herein, a computer-implemented method is provided comprising, using a processor of a device to solve for eigenvalues of a matrix system. The method comprises, using a processor of the computing device loading the matrix system into a memory of the computing device, and performing a domain decomposition into non-overlapping subdomains and a reordering of matrices of the matrix system. An interface variable projection subspace associated with interface variables of an adjacency graph of the matrix system is created, wherein the interface variables are related to nodes of the adjacency graph which are connected with nodes located in neighboring partitions. The method further comprises creating an internal variable projection subspace associated with internal variables of the adjacency graph of the matrix system, wherein the internal variables are related to nodes of the adjacency graph which are connected only to nodes located in the same partition. The method further comprises building a projection matrix based on the interface variable projection subspace and the internal variable projection subspace, determining eigenvalues that solve a Raleigh-Ritz eigenvalue problem utilizing the projection matrix, and utilizing the determined eigenvalues in a technical application system.

Advantageously, this approach may: a) preserve the advantages of algebraic domain decomposition eigenvalue solvers, such as reduced orthogonalization costs and inherent parallelism, while, at the same time, b) increase their accuracy without considering more than one shift.

According to another aspect disclosed herein, method may be applied to a technical application system that is selected from the group consisting of a clustering analytic system and a frequency response analytic system. Advantageously, when applying the techniques in these situations, a reduction in time and processing resources may be realized.

According to another aspect disclosed herein, the method may incorporate the clustering analytic system in an artificial intelligence (AI) classification system, and the frequency response analytic system may be selected from the group consisting of an audio analytic system and a building structure analytic system. The application to these technological areas may help reduce costs by requiring less computational demand.

According to another aspect disclosed herein, the method may further, when creating an interface variable projection subspace associated with interface variables, utilize a zeroth-order truncation of a spectral Schur complement. Advantageously, this avoids computation of a derivative of the spectral Schur complement.

According to another aspect disclosed herein, an apparatus for a technical system comprises a memory and a processor. The processor is configured for loading the matrix system into a memory of the computing device, and performing a domain decomposition into non-overlapping subdomains and a reordering of matrices of the matrix system. An interface variable projection subspace associated with interface variables of an adjacency graph of the matrix system is created, wherein the interface variables are related to nodes of the adjacency graph which are connected with nodes located in neighboring partitions. The processor is further configured for creating an internal variable projection subspace associated with internal variables of the adjacency graph of the matrix system, wherein the internal variables are related to nodes of the adjacency graph which are connected only to nodes located in the same partition. The processor is further configured for building a projection matrix based on the interface variable projection subspace and the internal variable projection subspace, determining eigenvalues that solve a Raleigh-Ritz eigenvalue problem utilizing the projection matrix, and utilizing the determined eigenvalues in a technical application system.

Advantageously, this approach may: a) preserve the advantages of algebraic domain decomposition eigenvalue solvers, such as reduced orthogonalization costs and inherent parallelism, while, at the same time, b) increase their accuracy without considering more than one shift.

According to another aspect disclosed herein, the apparatus may further, when creating an interface variable projection subspace associated with interface variables, utilize a zeroth-order truncation of a spectral Schur complement. Advantageously, this avoids computation of a derivative of the spectral Schur complement.

According to another aspect disclosed herein, the apparatus may be applied to a technical application system that is selected from the group consisting of a clustering analytic system and a frequency response analytic system. Advantageously, when applying the techniques in these situations, a reduction in time and processing resources may be realized.

According to another aspect disclosed herein, the apparatus may incorporate the clustering analytic system in an artificial intelligence (AI) classification system, and the frequency response analytic system may be selected from the group consisting of an audio analytic system and a building structure analytic system. The application to these technological areas may help reduce costs by requiring less computational demand.

According to another aspect disclosed herein, a computer program product is provided for a technical analytical system. The computer program product comprises one or more computer readable storage media, and program instructions collectively stored on the one or more computer readable storage media. The program instructions comprise program instructions for loading the matrix system into a memory of the computing device, and performing a domain decomposition into non-overlapping subdomains and a reordering of matrices of the matrix system. An interface variable projection subspace associated with interface variables of an adjacency graph of the matrix system is created, wherein the interface variables are related to nodes of the adjacency graph which are connected with nodes located in neighboring partitions. The instructions further configure the processor for creating an internal variable projection subspace associated with internal variables of the adjacency graph of the matrix system, wherein the internal variables are related to nodes of the adjacency graph which are connected only to nodes located in the same partition. The instructions further configure the processor for building a projection matrix based on the interface variable projection subspace and the internal variable projection subspace, determining eigenvalues that solve a Raleigh-Ritz eigenvalue problem utilizing the projection matrix, and utilizing the determined eigenvalues in a technical application system.

Advantageously, this approach may: a) preserve the advantages of algebraic domain decomposition eigenvalue solvers, such as reduced orthogonalization costs and inherent parallelism, while, at the same time, b) increase their accuracy without considering more than one shift.

DETAILED DESCRIPTION

Overview

The following acronyms may be used below:API application program interfaceARM advanced RISC machineCD-ROM compact disc ROMCMS content management systemCoD capacity on demandCPU central processing unitCUoD capacity upgrade on demandDPS data processing systemDVD digital versatile diskEPROM erasable programmable read-only memoryFPGA field-programmable gate arraysHA high availabilityIaaS infrastructure as a serviceI/O input/outputIPL initial program loadISP Internet service providerISA instruction-set-architectureLAN local-area networkLPAR logical partitionPaaS platform as a servicePDA personal digital assistantPLA programmable logic arraysRAM random access memoryRISC reduced instruction set computerROM read-only memorySaaS software as a serviceSLA service level agreementSRAM static random-access memoryWAN wide-area network
Data Processing System in General

FIG.1Ais a block diagram of an example DPS according to one or more embodiments. In this illustrative example, the DPS10may include communications bus12, which may provide communications between a processor unit14, a memory16, persistent storage18, a communications unit20, an I/O unit22, and a display24.

The processor unit14serves to execute instructions for software that may be loaded into the memory16. The processor unit14may be a number of processors, a multi-core processor, or some other type of processor, depending on the particular implementation. A number, as used herein with reference to an item, means one or more items. Further, the processor unit14may be implemented using a number of heterogeneous processor systems in which a main processor is present with secondary processors on a single chip. As another illustrative example, the processor unit14may be a symmetric multi-processor system containing multiple processors of the same type.

The memory16and persistent storage18are examples of storage devices26. A storage device may be any piece of hardware that is capable of storing information, such as, for example without limitation, data, program code in functional form, and/or other suitable information either on a temporary basis and/or a permanent basis. The memory16, in these examples, may be, for example, a random access memory or any other suitable volatile or non-volatile storage device. The persistent storage18may take various forms depending on the particular implementation.

For example, the persistent storage18may contain one or more components or devices. For example, the persistent storage18may be a hard drive, a flash memory, a rewritable optical disk, a rewritable magnetic tape, or some combination of the above. The media used by the persistent storage18also may be removable. For example, a removable hard drive may be used for the persistent storage18.

The communications unit20in these examples may provide for communications with other DPSs or devices. In these examples, the communications unit20is a network interface card. The communications unit20may provide communications through the use of either or both physical and wireless communications links.

The input/output unit22may allow for input and output of data with other devices that may be connected to the DPS10. For example, the input/output unit22may provide a connection for user input through a keyboard, a mouse, and/or some other suitable input device. Further, the input/output unit22may send output to a printer. The display24may provide a mechanism to display information to a user.

Instructions for the operating system, applications and/or programs may be located in the storage devices26, which are in communication with the processor unit14through the communications bus12. In these illustrative examples, the instructions are in a functional form on the persistent storage18. These instructions may be loaded into the memory116for execution by the processor unit14. The processes of the different embodiments may be performed by the processor unit14using computer implemented instructions, which may be located in a memory, such as the memory16. These instructions are referred to as program code38(described below) computer usable program code, or computer readable program code that may be read and executed by a processor in the processor unit14. The program code in the different embodiments may be embodied on different physical or tangible computer readable media, such as the memory16or the persistent storage18.

The DPS10may further comprise an interface for a network29. The interface may include hardware, drivers, software, and the like to allow communications over wired and wireless networks29and may implement any number of communication protocols, including those, for example, at various levels of the Open Systems Interconnection (OSI) seven layer model.

FIG.1Afurther illustrates a computer program product30that may contain the program code38. The program code38may be located in a functional form on the computer readable media32that is selectively removable and may be loaded onto or transferred to the DPS10for execution by the processor unit14. The program code38and computer readable media32may form a computer program product30in these examples. In one example, the computer readable media32may be computer readable storage media34or computer readable signal media36. Computer readable storage media34may include, for example, an optical or magnetic disk that is inserted or placed into a drive or other device that is part of the persistent storage18for transfer onto a storage device, such as a hard drive, that is part of the persistent storage18. The computer readable storage media34also may take the form of a persistent storage, such as a hard drive, a thumb drive, or a flash memory, that is connected to the DPS10. In some instances, the computer readable storage media34may not be removable from the DPS10.

Alternatively, the program code38may be transferred to the DPS10using the computer readable signal media36. The computer readable signal media36may be, for example, a propagated data signal containing the program code38. For example, the computer readable signal media36may be an electromagnetic signal, an optical signal, and/or any other suitable type of signal. These signals may be transmitted over communications links, such as wireless communications links, optical fiber cable, coaxial cable, a wire, and/or any other suitable type of communications link. In other words, the communications link and/or the connection may be physical or wireless in the illustrative examples.

In some illustrative embodiments, the program code38may be downloaded over a network to the persistent storage18from another device or DPS through the computer readable signal media36for use within the DPS10. For instance, program code stored in a computer readable storage medium in a server DPS may be downloaded over a network from the server to the DPS10. The DPS providing the program code38may be a server computer, a client computer, or some other device capable of storing and transmitting the program code38.

The different components illustrated for the DPS10are not meant to provide architectural limitations to the manner in which different embodiments may be implemented. The different illustrative embodiments may be implemented in a DPS including components in addition to or in place of those illustrated for the DPS10.

Cloud Computing in General

Characteristics are as Follows

Service Models are as Follows

Deployment Models are as Follows

Referring now toFIG.1B, illustrative cloud computing environment52is depicted. As shown, cloud computing environment52includes one or more cloud computing nodes50with which local computing devices used by cloud consumers, such as, for example, personal digital assistant (PDA) or cellular telephone54A, desktop computer54B, laptop computer54C, and/or automobile computer system54N may communicate. Nodes50may communicate with one another. They may be grouped (not shown) physically or virtually, in one or more networks, such as Private, Community, Public, or Hybrid clouds as described hereinabove, or a combination thereof. This allows cloud computing environment52to offer infrastructure, platforms and/or software as services for which a cloud consumer does not need to maintain resources on a local computing device. It is understood that the types of computing devices54A-N shown inFIG.1Bare intended to be illustrative only and that computing nodes50and cloud computing environment52can communicate with any type of computerized device over any type of network and/or network addressable connection (e.g., using a web browser).

Referring now toFIG.1C, a set of functional abstraction layers provided by cloud computing environment52(FIG.1B) is shown. It should be understood in advance that the components, layers, and functions shown inFIG.1Care intended to be illustrative only and embodiments of the invention are not limited thereto. As depicted, the following layers and corresponding functions are provided:

Any of the nodes50in the computing environment52as well as the computing devices54A-N may be a DPS10.

Determination and Use of Eigenvalues and Eigenvectors of Large-Scale Systems

A technical analytic system is described herein that produces a solution that may be subsequently applied to solve real-word technical problems. Such a technical analytic system may, for example, be a spectral clustering analytic system, a frequency response analytic system, or any other system that utilizes matrices as described herein. More generally, the system described herein can be seen as a computationally efficient approach to generate spectral embeddings. These embeddings can be used for the purposes of dimensionality reduction.

Spectral Clustering Analysis

Spectral clustering is a technique with roots in graph theory, where the approach is used to identify communities of nodes in a graph based on the edges connecting them. Clustering is one of the main tasks in unsupervised machine learning. The goal is to assign unlabeled data to groups, with an aim that similar data points get assigned to the same group. Spectral clustering is a popular variant of clustering which uses information from the eigenvalues (spectrum) of special matrices built from a graph or data set. In particular, spectral clustering applies clustering algorithms to identify different communities in which the data lie, but it does so on an embedded or sketched version of the initial data. This step, which can be also seen as a dimensionality reduction idea, is fulfilled through the computation of a few of the eigenvectors of a similarity matrix.

In multivariate statistics and the clustering of data, spectral clustering techniques make use of the spectrum (eigenvalues) of a similarity matrix of the data to perform dimensionality reduction before clustering in fewer dimensions. The similarity matrix is provided as an input and comprises of a quantitative assessment of the relative similarity of each pair of points in the dataset. In application to image segmentation, spectral clustering is known as segmentation-based object categorization.

Given an enumerated set of data points, the similarity matrix may be defined as a symmetric matrix A, where Aij≥0 represents a measure of the similarity between data points with indices i and j. The general approach to spectral clustering is to use a standard clustering method (there are many such methods, k-means is discussed below) on relevant eigenvectors of a Laplacian matrix of A. There are many different ways to define a Laplacian which have different mathematical interpretations, and so the clustering will also have different interpretations. The eigenvectors that are relevant are the ones that correspond to smallest several eigenvalues of the Laplacian except for the smallest eigenvalue which will have a value of zero, if the graph is connected. If the graph is not connected, then the number of zero eigenvalues is equal to the number of connected components.

Spectral clustering may relate to partitioning of a mass-spring system, where each mass is associated with a data point and each spring stiffness corresponds to a weight of an edge describing a similarity of the two related data points. Specifically, the eigenvalue problem describing transversal vibration modes of a mass-spring system is exactly the same as the eigenvalue problem for the graph Laplacian matrix defined as
L:=D−A,where D is the diagonal matrix

The masses that are tightly connected by the springs in the mass-spring system move together from the equilibrium position in low-frequency vibration modes, so that the components of the eigenvectors corresponding to the smallest eigenvalues of the graph Laplacian can be used for meaningful clustering of the masses. One related spectral clustering technique is the normalized cuts algorithm or Shi-Malik algorithm commonly used for image segmentation. It partitions points into two sets (B1, B2) based on the eigenvector v corresponding to the second-smallest eigenvalue of the symmetric normalized Laplacian defined as
Lnorm:=I−D−1/2AD−1/2

A mathematically equivalent algorithm takes the eigenvector corresponding to the largest eigenvalue of the random walk normalized adjacency matrix P=D−1A.

Knowing the eigenvectors, partitioning may be done in various ways, such as by computing the median m of the components of the second smallest eigenvector v, and placing all points whose component in v is greater than m in B1and the rest in B2. The algorithm can be used for hierarchical clustering by repeatedly partitioning the subsets in this fashion.

Spectral clustering is computationally expensive unless a fast technique to compute a large number of eigenvectors (i.e., spectral embeddings) is available. State-of-the-art approaches rely on the application of the (unrestarted) Lanczos method. Their cost scales at best linearly with respect to the input size but with a multiplicative factor, which depends quadratically on the number of eigenvectors sought. This becomes a major bottleneck when spectral clustering needs to be applied to large-scale datasets.

Frequency Response Analysis

Frequency response analysis is the quantitative measure of the output spectrum of a system or device in response to a stimulus, and is used to characterize the dynamics of the system. It is a measure of magnitude and phase of the output as a function of frequency, in comparison to the input. In simplest terms, if a sine wave is injected into a system at a given frequency, a linear system will respond at that same frequency with a certain magnitude and a certain phase angle relative to the input. A sinusoidal test signal may be used to measure points on the frequency response of a transfer function or impedance function.

Estimating the frequency response for a physical system generally involves exciting the system with an input signal, measuring both input and output time histories, and comparing the two through a process such as the Fast Fourier Transform (FFT). The frequency content of the input signal must cover the frequency range of interest because the results will not be valid for the portion of the frequency range not covered.

The frequency response of a system can be measured by applying a test signal, for example:applying an impulse to the system and measuring its response (see impulse response);sweeping a constant-amplitude pure tone through the bandwidth of interest and measuring the output level and phase shift relative to the input;applying a signal with a wide frequency spectrum (for example multifrequency signals (nonorthogonal frequency-discrete multiplexing of signals (N-OFDM or as the same SEFDM) and OFDM), digitally-generated maximum length sequence noise, or analog filtered white noise equivalent, like pink noise), and calculating the impulse response by deconvolution of this input signal and the output signal of the system.

The frequency response is characterized by the magnitude of the system's response, typically measured in decibels (dB) or as a decimal, and the phase, measured in radians or degrees, versus frequency in radians/sec or Hertz (Hz). These response measurements can be plotted in three ways: by plotting the magnitude and phase measurements on two rectangular plots as functions of frequency to obtain a Bode plot; by plotting the magnitude and phase angle on a single polar plot with frequency as a parameter to obtain a Nyquist plot; or by plotting magnitude and phase on a single rectangular plot with frequency as a parameter to obtain a Nichols plot. For audio systems with nearly uniform time delay at all frequencies, the magnitude versus frequency portion of the Bode plot may be all that is of interest. For the design of control systems, any of the three types of plots (Bode, Nyquist, Nichols) can be used to infer closed-loop stability and stability margins (gain and phase margins) from the open-loop frequency response, provided that for the Bode analysis the phase-versus-frequency plot is included. The form of frequency response for digital systems (as example FFT filters) are periodical with multiple main lobes and sidelobes.

If the system under investigation is nonlinear, then applying purely linear frequency domain analysis will not reveal all the nonlinear characteristics. To overcome these limitations, generalized frequency response functions and nonlinear output frequency response functions have been defined that allow the user to analyze complex nonlinear dynamic effects. The nonlinear frequency response methods reveal complex resonance, inter modulation, and energy transfer effects that cannot be seen using a purely linear analysis and are becoming increasingly important in a nonlinear world.

One important field where frequency response analysis holds a prominent role is in structural mechanics where complex structures are modeled through Partial Differential Equations. Because of the high-cost of conventional approaches, partially stemming from large frequency intervals, it is imperative to use techniques which are more efficient in terms of memory consumption and wall-clock time requirements. One such approach is to perform a modal analysis instead, where the entire structure is projected (embedded) onto a lower-dimensional subspace capturing most of the energy of the problem. For example, an object might vibrate mostly along the modes corresponding to the few lowest frequencies, thus justifying the computation of only those eigenvectors associated with the lowest eigenvalues of the matrix pencil associated with a Finite Element discretization of the complex structure.

Computation of Eigenvectors

Both the spectral clustering analytic system and the frequency response analytic system (and other analytic systems) depend on determining and applying eigenvectors (outputs), of a large matrix pencil (input) that routinely appears in these types of analytic systems. Described herein is a system to determine and apply such eigenvalues/eigenvectors. The disclosed system allows for the practical computation and application of eigenvalues of matrices which may be so large they cannot fit in system memory.

FIG.2Ais a block diagram that illustrates the general concept of such a system200described above. The technical analytic system210described herein may be run on a DPS10, and may comprise a processor212and memory214that may be used by an eigenvalue/eigenvector determiner230to execute various functions described herein. The technical analytic system210may utilize, as its input, a large sparse symmetric matrix pencil (A, M)205that is generated e.g., in a spectral clustering analysis system or a frequency response analysis system, as described above. Solving for the eigenvalues, eigenvectors, or eigenpairs of the matrix system allows an application to a technical application system250, such as, e.g., one that uses a spectral clustering analysis as might be used, for example, in AI classification application, or in a technical application system250such as, e.g., one that uses a frequency response analysis as might be used, for example, response to an applied impulse signal, multifrequency signals (nonorthogonal frequency-discrete multiplexing of signals (N-OFDM or as the same SEFDM) and OFDM), impulse responses by deconvolution of the input signal and the output signal of the system, and audio analytic system, and a building structure analytic system.

The symmetric generalized eigenvalue problem is formally defined as:
Ax=λMx.

In general, for a matrix A, if there exists a vector x which is not all 0s and a scalar λ such that Ax=λMx, then x is said to be an eigenvector of A with corresponding eigenvalue λ. The matrix A may be thought of as a function that maps vectors to new vectors. Most vectors will end up somewhere completely different when A is applied to them, but eigenvectors only change in magnitude. If one draws a line through the origin and the eigenvector, then after the mapping, the eigenvector would still land on the line. The amount which the vector is scaled along the line depends on λ. Eigenvectors help describe the dynamics of systems represented by matrices, including spectral clustering and frequency response analysis.

Graphs are a natural way to represent many types of data. A graph is a set of nodes with a corresponding set of edges which connect the nodes. The edges may be directed or undirected and can even have weights associated with them. For example, a network of routers on the internet can easily be represented as a graph. The routers are the nodes, and the edges are the connections between pairs of routers. Some routers might only allow traffic in one direction, so the edges could be directed to represent which direction traffic can flow. The weights on the edges could represent the bandwidth available along that edge. With this setup, one could then query the graph to find efficient paths for transmitting data from one router to another across the network.

Since spectral clustering requires the computation of the eigenvectors associated with a few smallest eigenvalues of a graph Laplacian matrix, given the large size of modern data collections, the eigenvector computation task can be extremely costly in terms of resources, such as processor bandwidth and memory.

Additionally, because of the large size of some of the structures undergoing frequency response analysis, solving the equation by classical approaches is extremely costly. An alternative is the modal superposition method, which requires the computation of eigenvectors associated with a large number of smallest eigenvalues.

The present system may be used to determine eigenvector derivatives of Schur complements in a novel manner. Known techniques, more importantly, give less accurate results than the technique disclosed herein and also require more memory.

Matrix Pencil

In linear algebra, if A0, A1, . . . , Alare n×n complex matrices for some nonnegative integer I, and Al≠0 (the zero matrix), then the matrix pencil of degree I is the matrix-valued function defined on the complex numbers

A particular case is a linear matrix pencil A−λB with λ ∈(or) where A and B are complex (or real) n×n matrices. This is denoted briefly with the notation (A, B).

A pencil is called regular if there is at least one value of λ such that det(A−AB)≠0 Eigenvalues of a matrix pencil (A, B) are all complex numbers λ for which det(A−AB)=0 The set of the eigenvalues is called the spectrum of the pencil and is written σ(A, B). Moreover, the pencil is said to have one or more eigenvalues at infinity if B has one or more 0 eigenvalues.

The present disclosure concerns a Rayleigh-Ritz projection scheme to determine a selected number of eigenvalues (and, optionally, associated eigenvectors) of large and sparse symmetric pencils. This represents an improvement over the traditional shift-and-invert Lanczos method previously used. The focus lies in the determination of a few eigenvalues located immediately on the right (or left) of some real scalar, and the proposed technique is the more effective when a large number of eigenvalues/eigenvectors are sought. Such eigenvalue problems appear in applications such as low-frequency response analysis and spectral clustering, among others.

Determining eigenvalues located the closest to a given scalar is typically achieved by enhancing the projection method of choice by shift-and-invert. When the required accuracy in the approximation of the sought eigenvalues is not high, an alternative to shift-and-invert is the Automated Multilevel Substructuring (AMLS) technique. AMLS builds a large projection subspace while avoiding excessive orthogonalization costs and can be considerably faster than shift-and-invert Krylov subspace techniques in applications where a large number of eigenvalues are sought.

Process Flow

FIG.3is a flowchart that illustrates a process300for finding eigenpairs (eigenvalues, eigenvectors) according to some embodiments.

In operation305, the technical analytic system210may load/assemble the matrix system, e.g., the large sparse symmetric matrix pencil (A, M)205, which constitutes the domain. The present focus here assumes that matrices A and M are symmetric, and M is a symmetric positive definite (SPD) matrix. The matrix pencil (A, M)205has n eigenpairs, denoted here by (λi, x(i)), i=1, . . . , n, where λiconstitute eigenvalues, and x(i)constitute eigenvectors. Of interest is the determination of specified neveigenpairs for matrix pencil (A, M)205, which are located immediately on the right (or left) of some scalar α ∈.FIG.2Billustrates the set of real numbers R280showing the real axis282and the scalar α284, and the eigenvalues λi286.FIG.2Cillustrates a graph290comprising interior nodes292within partitions, and interface nodes294between partitions.

The loading/assembling of the matrix system includes inputting the matrices A and M, the scalars nev, which is the number of eigenpairs to determine, and α, which is a scalar value from which the eigenvalues are sought. It also includes inputting a number of partitions p, and a number of local modes computing from each substructure (parition).

The standard approach to compute a few eigenpairs of sparse and symmetric matrix pencils is by applying the Rayleigh-Ritz technique onto some carefully constructed subspaceofn. The goal is to find a subspacethat includes an invariant subspace associated with the nevsought eigenvalues λ1, . . . , λnev. The sought eigenpairs can then be identified (in the absence of roundoff errors) as a subset of the Ritz pairs of the matrix pencil (ZTAZ ZTMZ), where the matrix Z represents a basis of the subspace.

The matrices

A=(B1E1B2E2⋱⋮BPEPE1TE2T⋯EPTC),M=(MB(1)ME(1)MB(2)ME(2)⋱⋮MB(p)ME(p)(ME(1))T(ME(2))T⋯(ME(p))TMC)
may be written in a 2×2 block form:

A=(BEETC),M=(MBMEMETMC),
where B and MBare square matrices of size d×d, E and MEare rectangular matrices of size d×s, C and MCare square matrices of size s×s, n−d+s, and s≥nev. Similarly, the eigenvector x(i)associated with eigenvalue λimay be rewritten in a 2×1 block form

Ax(i)=λiMx(i)using the block form gives

(A-λi⁢M)⁢x(i)=(B-λi⁢MBE-λi⁢MEET-λi⁢METC-λi⁢MC)⁢(u(i)y(i))=0.
Eliminating u(i) from the second row equation leads to the s×s nonlinear eigenvalue problem

[C-λi⁢MC-(E-λi⁢ME)T⁢(B-λi⁢MB)-1⁢(E-λi⁢ME)︸S⁡(·)]⁢y(i)=0,u(i)=-(B-λi⁢MB)-1⁢(E-λi⁢ME)︸block-diagonal⁢y(i),
from which the eigenpair can be determined. The missing d×1 part of x(i), u(i)is then recovered by solving the linear system
(B−λiMB)u(i)=−(E−λiME)y(i).

From a domain decomposition perspective, an ideal choice is to set=u⊕y, and perform a Rayleigh-Ritz projection onto this subspace, where
y=span{y(i)}i=1, . . . ,nev,
u=span{−(B−λiMB)−1(E−λiME)y(1)}i=1, . . . ,nev,

The main goal of algebraic domain decomposition eigenvalue solvers is to build a projection subspace which, ideally, includes the subspace Z.

Domain Decomposition (Partitioning) into Subdomains and Reordering of Matrices

The first operation performed on the matrix pencil (A, M)205, referring toFIG.2D, is having the eigenvalue/eigenvector determiner230performing a domain decomposition (partitioning) into subdomains and reordering the matrices, in operation310. A graph associated with a given matrix pencil may be partitioned into a number of non-overlapping subdomains with the help of an algebraic graph partitioner232. The projection subspace may have interface variables234and interior variables236.

Practical applications of algebraic domain decomposition eigenvalue solvers rely on relabeling the unknowns/equations of the eigenvalue problem Ax=λMx such that matrices B−λMBand E−λMEare block-diagonal. This can be easily achieved by applying a graph partitioner to the adjacency graph of the matrix |A|+|M|. In this disclosure, only p-way partitionings are considered, and the partitioner divides the graph into p>1 non-overlapping subdomains. The rows/columns of matrices A and M are then reordered so that unknowns/equations associated with interior variables (i.e., nodes of the adjacency graph which are connected only to nodes located in the same partition) are listed before those associated with interface variables (i.e., nodes of the adjacency graph which are connected with nodes located in neighboring partitions). The sought eigenpairs of the original pencil are then approximated by a Rayleigh-Ritz projection on a subspace formed by the direct sum of the subspaces associated with the two different types of variables.

Algebraic Domain Decomposition

This disclosure focuses on algebraic domain decomposition eigenvalue solvers where the concept of domain decomposition is applied directly to the eigenvalue equation. The main challenge of algebraic domain decomposition eigenvalue solvers is the construction of the interface projection subspace due to the non-linear nature of the interface matrix operator, also known as “spectral Schur complement”.

In linear algebra and the theory of matrices, the Schur complement of a block matrix is defined as follows. Suppose A, B, C, D are respectively p×p, p×q, q×p, and q×q matrices, and D is invertible. Let

M=[ABCD]
so that M is a (p+q)×(p+q) matrix.

Then the Schur complement of the block D of the matrix M is the p×p matrix defined by
M/D:=A−BD−1C
and, if A is invertible, the Schur complement of the block A of the matrix M is the q×q matrix defined by
M/A:=C−CA−1B.

In the case that A or D is singular, substituting a generalized inverse for the inverses on M/A and M/D yields the generalized Schur complement.

AMLS, an alternate technique to that focused on here, can be seen as building this subspace by computing a few of the eigenvectors of the linear generalized eigenvalue problem stemming by a first-order truncation of the Taylor series of the spectral Schur complement expanded around some real shift. As a result, the accuracy provided by AMLS might deteriorate considerably as one moves further away from the chosen shift. An alternative suggested is to form the interface projection subspaces by considering a complex rational transformation to the original pencil which damps the interface components of those eigenvectors associated with unwanted eigenvalues. While these techniques can indeed lead to enhanced accuracy, they also require multiple matrix factorizations computed in complex arithmetic, which is processor intensive.

The algorithm disclosed herein may: a) preserve the advantages of algebraic domain decomposition eigenvalue solvers, such as reduced orthogonalization costs and inherent parallelism, while, at the same time, b) increase their accuracy without considering more than one shift. Numerical experiments have demonstrated the efficiency of the proposed technique on sequential/distributed memory architectures as well as its competitiveness against schemes such as shift-and-invert Lanczos, and AMLS combined with p-way vertex-based partitionings.

Method Characteristics

Primary characteristics of the method include: 1) zeroth-order truncation of the interface matrix operator; 2) exploiting Taylor series expansions of interface eigenvectors; and 3) reduced orthogonalization costs and enhanced parallelism. These are summarized below.

Zeroth-Order Truncation of the Interface Matrix Operator

In contrast to AMLS, the algorithm disclosed herein only considers a zeroth-order truncation of the spectral Schur complement (the interface matrix operator), i.e., there is only a need to partially solve a standard eigenvalue problem to generate the interface projection subspace234. This approach avoids the need to compute/apply the derivative of the spectral Schur complement, as AMLS does.

Exploiting Taylor Series Expansions of Interface Eigenvectors

The accuracy of the interface projection subspace is enhanced by expanding the (analytic) eigenvectors of the spectral Schur complement through their Taylor series, and injecting a few leading eigenvector derivatives into the subspace. It is shown theoretically, and verified experimentally, that injecting up to second leading eigenvector derivatives leads to eigenvalue approximation errors of up to quartic order. These eigenvector derivatives are computed independently of each other by exploiting deflated Krylov subspace solvers. Several details are discussed, and bounds on the effective condition numbers are given.

Reduced Orthogonalization Costs and Enhanced Parallelism

Similarly to domain decomposition schemes such as AMLS and RF-DDES, orthonormalization is only applied to vectors whose length is equal to the number of interface variables instead of the entire set of domain variables (e.g., as in shift-and-invert Lanczos). This becomes important when both a large number of eigenvalues is sought and the number of interface variables is much smaller compared to the total number of equations/unknowns.

In addition, the projection subspace associated with the interior variables of each subdomain is built independently of each other (and thus trivially parallel) by computing local eigenmodes and computing up to second-order resolvent expansions. Details on the expected accuracy of the eigenvalue approximations are presented.

This approach also has the advantages of parallelism, and orthogonalization costs, and it is derivative-free, and may yield greater accuracy, although it may or may not reduce computational cost.

The computation of dϕyψi(α) requires a linear system solution with matrix S(α)−μψi(α)|. It is possible to take advantage of the computed eigenvectors of S(α) (deflation). Computing more than one or two leading derivatives is rather impractical.

FIG.2Eis a graph240that illustrates the elements within suboperation (1) of operation315of utilizing truncated spectral Schur complements. Let μj(ζ) denote the jtheigenvalue of the matrix. Then:

S⁡(ζ)=C-ζ⁢MC-∑i=1d(E-ζ⁢ME)T⁢vi⁢vjT(E-ζ⁢ME)δi-ζ,
where (δi) denotes the itheigenpair of (B, MB). The eigenvalues δi∈ Λ (B, MB) are poles of S(ζ), and the eigenvalue curves are analytic except at their poles.

Assuming now that μj(A)=0 where j is known, and let σ ∈such that σ ∈ [λ−, λ+] includes only the root of λ of μj(.). Ideally, μj(σ)=0, i.e., σ ≡ λ and yj(λ) is an eigenvector of S(ζ)y=0. However, in practice, it may or may not be true that σ≈λ, so the focus is on how to improve upon the span {yj(σ)}.

Expanding yj(σ) through its Taylor series around σ gives:

yj(λ)=∑i=0∞(λ-σ)ii!⁢(di⁢yj(ζ)d⁢ζi)ζ=σ.
The eigenvector x associated with λ may be written as

Referring to the graph240′ inFIG.2F, the above discussion omits several practical details, including: a) it assumes A is known; b) it assumes the subscript j for which μj(λ)=0 is known; and c) it requires formulas to approximate y, dy, d2y, . . . . In practice, it is expected that consecutive eigenvalue curves will cross the real axis at consecutive λ∈Λ (A, M). Thus the algorithm needs to compute only the neveigenvalue curves associated with the nevsmallest non-negative eigenvalues of S(α).

Operation320creates a projection of the subspace associated with the interior variables. The main idea here is to write Bλ−1=Bσ−1Σi=0p[(λ−σ)MBBσ−1]iand retain a first-order approximation plus eigenvectors located the closest to σ. Combining p≥1 with μ=1, it can be shown that
|λ−{circumflex over (λ)}|=O((λ−σ)χ), where
χ=min{4.2(τ+1)·τ+3}.

The operation takes as input, σ:=α, κ, and Y:=orthonormal basis returned by the previous (interface variables) algorithm. It then loops for j=1, . . . , p:
Vj=eigs(Bα(j),MB(j),κ,sm)(Bα(j)=Bj−αMB(j)).
It then forms:

Next, in operation325, the projection matrix Z is built for the Rayleigh-Ritz eigenvalue problem as follows. The matrix Y returned from operation315is distributed among the p subdomains, and can be written as
Y=[Y1TY2T. . . YpT]T
where Yj∈sj×ηis the row block of matrix Y associated with the jthsubdomain, and η ∈* denotes the column dimension of matrix Y. By definition, η is equal to an integer multiple of nev. The following matrices are defined:

The projection matrix z can then be written as:

The total memory overhead associated with the jthsubdomain is equal to that of storing κdi+(3di+si)η floating-point scalars. The dimension of the Rayleigh-Ritz pencil (ZTAZ, ZTMZ) is equal to κP+3η and can be solved by the appropriate routine in LAPACK, e.g., dsygv. When ME(j)=0, Pj=Eα(j)Yj, and the bottom row-block of matrix Z becomes [0s,pκ, Y, 0s,η]. The dimension of that Rayleight-Ritz pencil then reduces to κp+2η.

Finally, in operation330, the Rayleigh-Ritz eigenvalue problem is solved by solving:
(ZTAZ){tilde over (x)}={tilde over (λ)}(ZTMZ){tilde over (x)}
for the nevsought ({tilde over (λ)}, {tilde over (x)}).

The following compares the presently disclosed algorithm with the shift-and-invert Lanczos when the latter is applied directly to the pencil (A—αM, M).

The first difference between these two techniques is orthogonalization cost. Applying k steps of shift-and-invert Lanczos to matrix pencils (S(α), I) and (A—aM, M) leads to a total orthogonalization cost of O(k2s) and O(k2n), respectively. Thus, the present algorithm reduces orthogonalization costs by a factor of n/s, and this difference becomes more pronounced as nevincreases (since k≥nev). The second difference between the present algorithm and shift-and-invert Lanczos is the number of linear system solutions with matrix Bα. It is straightforward to verify that applying k steps of shift-and-invert Lanczos to the pencil (A—aM, M) requires 2k such linear system solutions. In contrast, the present algorithm requires 3η (2η if ME=0) linear solves with Bα. Nonetheless, those 3η linear solves with matrix Bαmay be performed simultaneously, since all right-hand sides are available at the same time. Thus, the associated cost can be marginally higher than that of solving for a few right-hand sides.

To the above cost one also adds the computational cost required to compute the eigenvectors of the block-diagonal pencil (Bα, MB) in the present algorithm. On the other hand, the accuracy achieved by shift-and-invert Lanczos applied directly to pencil (A—aM, M) can be significantly higher than that of the present algorithm. Nonetheless, in many applications, the sought eigenpairs are not required up to a very high accuracy, e.g., eigenvalue problems originating from discretizations need be solved up to the associated discretization error.

The present algorithm is based on domain decomposition and is thus well-suited for execution in modern distributed memory environments. For example, each one of the p subdomains can be mapped to a separate message passing interface (MPI) process.

Performing any type of operations with the block-diagonal matrices

Bσ(j), Eσ(j), and ME(j)is then an entirely local process in the jthsubdomain (i.e., creating the projection of the subspace associated with the interior variables in operation320is extremely parallel). An additional layer of parallelism is possible in each subdomain by further exploiting shared memory parallelism to perform the required computations with the aforementioned matrices.

In contrast to operation320, operation312involves computations with the Schur complement matrix S(α), which is distributed among the p subdomains. In this case, point-to-point communication among neighboring subdomains is necessary. Finally, the Rayleigh-Ritz eigenvalue problem is typically small enough so that it can be replicated in all MPI processes and solved redundantly.

A detailed application of the eigenvectors of the pencil (A,M) computed by the present system on spectral clustering is as follows. Let A denote the graph Laplacian matrix of an undirected graph and M be a diagonal matrix where its ithdiagonal entry holds the degree of the ithvertex of the graph. Now let x(1), x(2), . . . , x(k), denote the eigenvectors associated with the k algebraically smallest eigenvalues of the pencil (A,M), and let the columns of matrix X be equal to the k computed eigenvectors x(1), x(2), . . . , x(k). Now let xidenote the ithrow of matrix X. The final step of spectral clustering is to call a clustering algorithm (e.g., k-means) to partition the k-dimensional points determined by the rows of matrix X into a user-given number of clusters.

Technical Application

The one or more embodiments disclosed herein accordingly provide an improvement to computer technology. For example, an improvement to an analytical computer that is able to efficiently solve the eigenvalue/eigenvector problem for large and sparce symmetric matrix pencils provides a technical advantage in terms of computer resources, such as processor cycles and memory.

Computer Readable Media

The present invention may be a system, a method, and/or a computer readable media at any possible technical detail level of integration. The computer program product may include a computer readable storage medium (or media) having computer readable program instructions thereon for causing a processor to carry out aspects of the present invention.