Determining and using spectral embeddings of large systems by substructuring.

Abstract

A device solves eigenvalues ​​of a matrix system. The device performs domain partitioning of the matrix system into non-overlapping subdomains and reordering of the matrices of the matrix system. A boundary variable projection subspace associated with boundary variables of an adjacency graph of the matrix system is created. The boundary variables are associated with nodes of the adjacency graph that are connected to nodes located in neighboring partitions. An internal variable projection subspace associated with internal variables of the adjacency graph of the matrix system is created, where the internal variables are associated with nodes of the adjacency graph that are connected only to nodes located in the same partition. A projection matrix is ​​constructed based on the boundary variable projection subspace and the internal variable projection subspace. The device uses the projection matrix to determine eigenvalues ​​that solve the Rayleigh-Ritz eigenvalue problem.

Description

[Background technology]

[0001] Disclosed herein are systems and associated methods for determining and utilizing eigenvectors (embeddings) of large sparse symmetric matrix pencils that may arise in a variety of technical applications including, for example, spectral clustering and frequency response analysis.

[0002] The present 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 a large sparse symmetric matrix pencil. These eigenpairs are computed via Rayleigh-Ritz projection. The projected subspace is divided into two separate subspaces, each one associated with one of two types of variables: interior and boundary. The projected subspace associated with the boundary variables is constructed by computing a small number of eigenvectors and associated first and second derivatives of a zeroth order approximation of a nonlinear matrix-valued boundary operator.

[0003] Projected subspaces associated with the internal variables are constructed independently in each subregion by exploiting local eigenmodes and matrix resolvent approximations. The sought eigenpairs are then approximated via Rayleigh-Ritz projection. The proposed technique may result in computational resource improvements over other processes such as shift-and-reverse Lanczos and automated multilevel substructuring combined with p-vertex-based partitioning. Summary of the Invention

[0004] According to one aspect disclosed herein, a computer-implemented method is provided that includes solving eigenvalues ​​of a matrix system using a processor of a device. The method includes using a processor of the computing device, loading the matrix system into a memory of the computing device, and performing a domain partition into non-overlapping subdomains and a reordering of matrices of the matrix system. A boundary variable projection subspace associated with boundary variables of an adjacency graph of the matrix system is created, where the boundary variables are related to nodes of the adjacency graph that are connected to nodes located in neighboring partitions. The method further includes creating an internal variable projection subspace associated with internal variables of the adjacency graph of the matrix system, where the internal variables are related to nodes of the adjacency graph that are only connected to nodes located in the same partition. The method further includes constructing a projection matrix based on the boundary variable projection subspace and the internal variable projection subspace, utilizing the projection matrix to determine eigenvalues ​​that solve a Rayleigh-Ritz eigenvalue problem, and utilizing the determined eigenvalues ​​in a technical application system.

[0005] Advantageously, this approach: a) preserves the advantages of algebraic domain decomposition eigensolvers, such as reduced orthogonalization cost and inherent parallelism, while at the same time b) may increase their accuracy without considering more than one shift.

[0006] According to another aspect disclosed herein, the methods may be applied to technical application systems selected from the group consisting of clustering analysis systems and frequency response analysis systems. Advantageously, savings in time and processing resources may be realized when applying these techniques in these contexts.

[0007] According to another aspect disclosed herein, the method may incorporate a clustering analysis system into an artificial intelligence (AI) classification system, and the frequency response analysis system may be selected from the group consisting of an audio analysis system and a building structural analysis system. Application to these technology areas may help reduce costs by requiring less computational demands.

[0008] According to another aspect disclosed herein, the method may further utilize a zeroth order truncation of the spectral Schur complement when creating the boundary variable projected subspace associated with the boundary variables, advantageously avoiding the calculation of derivatives of said spectral Schur complement.

[0009] According to another aspect disclosed herein, an apparatus for a technical system comprises a memory and a processor. The processor is configured to load a matrix system into the memory of the computing device and to perform a domain partitioning into non-overlapping subdomains and a reordering of matrices of the matrix system. A boundary variable projection subspace associated with boundary variables of an adjacency graph of the matrix system is created, where the boundary variables are related to nodes of the adjacency graph that are connected to nodes located in neighboring partitions. The processor is further configured to create an internal variable projection subspace associated with internal variables of the adjacency graph of the matrix system, where the internal variables are related to nodes of the adjacency graph that are only connected to nodes located in the same partition. The processor is further configured to construct a projection matrix based on the boundary variable projection subspace and the internal variable projection subspace, to utilize the projection matrix to determine eigenvalues ​​that solve a Rayleigh-Ritz eigenvalue problem, and to utilize the determined eigenvalues ​​in a technical application system.

[0010] Advantageously, this approach: a) preserves the advantages of algebraic domain decomposition eigensolvers, such as reduced orthogonalization cost and inherent parallelism, while at the same time b) may increase their accuracy without considering more than one shift.

[0011] According to another aspect disclosed herein, the apparatus may further utilize a zeroth order truncation of the spectral Schur complement when creating the boundary variable projected subspace associated with the boundary variables, advantageously avoiding the calculation of derivatives of said spectral Schur complement.

[0012] According to another aspect disclosed herein, the apparatus may be applied to said technical application systems selected from the group consisting of clustering analysis systems and frequency response analysis systems. Advantageously, savings in time and processing resources may be realized when applying these techniques in these contexts.

[0013] According to another aspect disclosed herein, the device may incorporate a clustering analysis system into an artificial intelligence (AI) classification system, and the frequency response analysis system may be selected from the group consisting of an audio analysis system and a building structure analysis system. Application to these technology areas may help reduce costs by requiring less computational demands.

[0014] According to another aspect disclosed herein, a computer program product for an engineering analysis system is provided. 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 include program instructions for loading a matrix system into a memory of the computing device, and program instructions for performing a domain partition into non-overlapping subdomains and a reordering of matrices of the matrix system. A boundary variable projected subspace associated with boundary variables of an adjacency graph of the matrix system is created, where the boundary variables are related to nodes of the adjacency graph that are connected to nodes located in neighboring partitions. The instructions further configure the processor to create an internal variable projected subspace associated with internal variables of the adjacency graph of the matrix system, where the internal variables are related to nodes of the adjacency graph that are connected only to nodes located in the same partition. The instructions further configure the processor to construct a projection matrix based on the boundary variable projected subspace and the internal variable projected subspace, to determine eigenvalues ​​that solve a Rayleigh-Ritz eigenvalue problem using the projection matrix, and to utilize the determined eigenvalues ​​in a technical application system.

[0015] Advantageously, this approach: a) preserves the advantages of algebraic domain decomposition eigensolvers, such as reduced orthogonalization cost and inherent parallelism, while at the same time b) may increase their accuracy without considering more than one shift.

[0016] Furthermore, embodiments may take the form of an associated computer program product accessible from a computer usable or computer readable medium providing program code for use by or in connection with a computer or any instruction execution system. For purposes of this description, a computer usable or computer readable medium may be any apparatus that may include a mechanism for storing, communicating, propagating or transferring a program for use by or in connection with an instruction execution system, apparatus, or device. [Brief description of the drawings]

[0017] Various embodiments are described herein with reference to different subject matter. In particular, some embodiments may be described with reference to methods, while other embodiments may be described with reference to devices and systems. However, a person skilled in the art will infer from the above and following description that, unless otherwise informed, any combination of features belonging to one type of subject matter, as well as any combination between features relating to different subject matters, in particular between method features and device and system features, is considered to be disclosed within this document.

[0018] The above defined aspects, and further aspects disclosed herein, will be apparent from and will be explained with reference to one or more example embodiments described hereinafter, to which the invention is not limited. Various embodiments are described, by way of example only, with reference to the following drawings, in which:

[0019] [Figure 1A] FIG. 1 is a block diagram of a data processing system (DPS) according to one or more embodiments disclosed herein.

[0020] [Figure 1B] FIG. 1 is a pictorial diagram illustrating a cloud computing environment according to one embodiment disclosed herein.

[0021] [Figure 1C] FIG. 2 is a pictorial diagram illustrating abstraction model layers according to one embodiment disclosed herein.

[0022] [Figure 2A] FIG. 1 is a block diagram illustrating a technical analysis system according to some implementations.

[0023] [Figure 2B] 11 is a graph showing eigenvalue solutions adjacent to a sought value of α according to some implementations.

[0024] [Figure 2C] 1 is a pictorial diagram showing interior nodes versus boundary nodes in an adjacency graph according to some implementations.

[0025] [Figure 2D] FIG. 2 is a block diagram illustrating an eigenvalue / eigenvector determiner according to some implementations.

[0026] [Figure 2E] 13 is a graph showing eigenvalue curves in subregions of those domains. [Figure 2F] 13 is a graph showing eigenvalue curves in subregions of those domains.

[0027] [Diagram 3] 1 is a flowchart illustrating a process used herein according to some implementations. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0028] Summary The following acronyms may be used below: API Application Programming Interface ARM Advanced RISC Machine CD-ROM Compact Disc ROM CMS Content Management System CoD Capacity on Demand CPU Central Processing Unit CUoD Capacity Upgrade On Demand DPS Data Processing System DVD Digital Versatile Disc EPROM Erasable Programmable Read Only Memory FPGA Field Programmable Gate Array HA high availability IaaS Infrastructure as a Service I / O Input / Output IPL Initial Program Load ISP Internet Service Provider ISA Instruction Set Architecture LAN Local Area Network LPAR Logical Partition PaaS Platform as a Service PDA Personal Digital Assistant PLA Programmable Logic Array RAM Random Access Memory RISC Reduced Instruction Set Computer ROM Read Only Memory SaaS Software as a Service SLA Service Level Agreement SRAM Static Random Access Memory WAN Wide Area Network General Data Processing System

[0029] 1A is a block diagram of an example DPS according to one or more embodiments. In this illustrative example, DPS 10 may include a communication bus 12 that may provide communication between a processor unit 14, a memory 16, persistent storage 18, a communication unit 20, an I / O unit 22, and a display 24.

[0030] Processor unit 14 functions to execute instructions for software that may be loaded into memory 16. Processor unit 14 may be multiple processors, a multi-core processor, or some other type of processor, depending on the particular implementation. When used herein with reference to an item, multiple means one or more items. Additionally, processor unit 14 may be implemented using multiple heterogeneous processor systems where a main processor resides on a single chip with secondary processors. As another illustrative example, processor unit 14 may be a symmetric multi-processor system that includes multiple processors of the same type.

[0031] Memory 16 and persistent storage 18 are examples of storage devices 26. A storage device may be any piece of hardware capable of storing information, such as, for example, without limitation, data, program code in a functional form, and / or other suitable information, on either a temporary and / or persistent basis. Memory 16, in these examples, may be, for example, a random access memory or any other suitable volatile or non-volatile storage device. Persistent storage 18 may take various forms depending on the particular implementation.

[0032] For example, persistent storage 18 may include one or more components or devices. For example, persistent storage 18 may 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 persistent storage 18 may be removable. For example, a removable hard drive may be used for persistent storage 18.

[0033] The communications unit 20 in these examples may provide for communications with other DPSs or devices. In these examples, the communications unit 20 is a network interface card. The communications unit 20 may provide for communications through the use of one or both of physical and wireless communications links.

[0034] Input / output unit 22 may allow for the input and output of data with other devices that may be connected to DPS 10. For example, input / output unit 22 may provide a connection for user input through a keyboard, mouse, and / or some other suitable input device. Additionally, input / output unit 22 may send output to a printer. Display 24 may provide a mechanism for displaying information to a user.

[0035] Instructions for the operating system, applications, and / or programs may be located in storage devices 26, which are in communication with processor unit 14 via communications bus 12. In these illustrative examples, the instructions are in a functional form on persistent storage 18. These instructions may be loaded into memory 116 for execution by processor unit 14. The processes of the different embodiments may be executed by processor unit 14 using computer-implemented instructions, which may be located in a memory, such as memory 16. These instructions are referred to as program code 38 (described below), computer usable program code, or computer readable program code, which may be read and executed by a processor in processor unit 14. The program code in the different embodiments may be embodied on different physical or tangible computer readable media, such as memory 16 or persistent storage 18.

[0036] DPS 10 may further include an interface for network 29. The interface may include hardware, drivers, software, etc. to enable communication over wired and wireless network 29 and may implement any number of communication protocols including, for example, protocols at various levels of the Open Systems Interconnection (OSI) seven-layer model.

[0037] 1A further illustrates computer program product 30, which may include program code 38. Program code 38 may be located in a functional form on a computer readable medium 32, which is selectively removable, and may be loaded onto or transferred to DPS 10 for execution by processor unit 14. Program code 38 and computer readable medium 32 may form computer program product 30, in these examples. In one example, computer readable medium 32 may be computer readable storage medium 34 or computer readable signal medium 36. Computer readable storage medium 34 may include, for example, an optical or magnetic disk that is inserted or placed into a drive or other device that is part of persistent storage 18 for transfer to a storage device, such as a hard drive that is part of persistent storage 18. Computer readable storage medium 34 may also take the form of persistent storage, such as a hard drive, thumb drive, or flash memory, that is connected to DPS 10. In some cases, computer readable storage medium 34 may not be removable from DPS 10.

[0038] Alternatively, program code 38 may be transferred to DPS 10 using computer readable signal medium 36. Computer readable signal medium 36 may be, for example, a propagated data signal containing program code 38. For example, computer readable signal medium 36 may 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, fiber optic cable, coaxial cable, a wire, and / or any other suitable type of communications link. In other words, in the illustrative examples, the communications links and / or connections may be physical or wireless.

[0039] In some exemplary embodiments, program code 38 may be downloaded over a network from another device or DPS to persistent storage 18 via computer readable signal medium 36 for use within DPS 10. For example, program code stored on a computer readable storage medium in a server DPS may be downloaded over a network from the server to DPS 10. The DPS providing program code 38 may be a server computer, a client computer, or some other device capable of storing and transmitting program code 38.

[0040] The different components illustrated for DPS 10 are not intended to provide architectural limitations to the manner in which different embodiments may be implemented. Different illustrative embodiments may be implemented in a DPS that includes components in addition to, or in place of, the components illustrated for DPS 10. General Cloud Computing

[0041] Although this disclosure includes detailed descriptions of cloud computing, it should be understood that implementation of the teachings described herein is not limited to a cloud computing environment. Rather, embodiments of the invention can be implemented in conjunction with any other type of computing environment, now known or later developed.

[0042] Cloud computing is a model of service delivery for enabling convenient, on-demand network access to a shared pool of configurable computing resources (e.g., networks, network bandwidth, servers, processing, memory, storage, applications, virtual machines, and services) that can be rapidly provisioned and released with minimal administrative effort or interaction with a service provider. The cloud model may include at least five characteristics, at least three service models, and at least four deployment models. The characteristics are as follows:

[0043] On-Demand Self-Service: Cloud consumers can unilaterally provision computing capacity, such as server time and network storage, automatically as needed, without the need for human interaction with the service provider.

[0044] Broad network access: The capabilities are available across the network and accessed through standard mechanisms that facilitate use by heterogeneous thin or thick client platforms (eg, cell phones, laptops, and PDAs).

[0045] Resource Pooling: Provider computing resources are pooled to serve multiple consumers using a multi-tenant model, with different physical and virtual resources dynamically allocated and reallocated according to demand. There is location independence in that consumers generally have no control or knowledge over the exact location of the provided resources, but it may be possible to specify the location at a higher level of abstraction (e.g., country, state, or data center).

[0046] Rapid Elasticity: This capacity can be rapidly and elastically provisioned, in some cases automatically, to rapidly scale out, and rapidly released to rapidly scale in. To the consumer, the capacity available for provisioning often appears unlimited, and can be purchased in any amount at any point in time.

[0047] Measured Services: Cloud systems automatically control and optimize resource usage by leveraging metering capabilities at a level of abstraction appropriate to the type of service (e.g., storage, processing, bandwidth, and active user accounts). Resource usage can be monitored, controlled, and reported, providing transparency to both providers and consumers of the services utilized. The service model is as follows:

[0048] Software as a Service (SaaS): The ability offered to the consumer is to use a provider's applications running on a cloud infrastructure. The applications are accessible from a variety of client devices through thin-client interfaces such as web browsers (e.g., web-based email). The consumer does not manage or control the underlying cloud infrastructure, including networks, servers, operating systems, storage, or even individual application capabilities, with the possible exception of limited user-specific application configuration settings.

[0049] Platform as a Service (PaaS): The ability offered to a consumer is to deploy applications that the consumer has created or acquired, written using programming languages ​​and tools supported by the provider, onto a cloud infrastructure. The consumer does not manage or control the underlying cloud infrastructure, including networks, servers, operating systems, or storage, but does control the deployed applications and, in some cases, the application hosting environment configuration.

[0050] Infrastructure as a Service (IaaS): The ability offered to consumers is to provision processing, storage, network and other basic computing resources on which the consumer can deploy and run any software, which may include operating systems and applications. The consumer does not manage or control the underlying cloud infrastructure, but does control the operating systems, storage, deployed applications, and possibly limited control over selected networking components (e.g., host firewalls). The deployment model is as follows:

[0051] Private Cloud: This cloud infrastructure operates solely for an organization. It may be managed by that organization or a third party and may exist on-premise or off-premise.

[0052] Community Cloud: This cloud infrastructure is shared by several organizations and supports a particular community with shared concerns (e.g., mission, security requirements, policies, and compliance considerations). The community cloud may be managed by those organizations or a third party and may exist on-premise or off-premise.

[0053] Public Cloud: The cloud infrastructure is made available to the general public or large industry entities and is owned by organizations that sell cloud services.

[0054] Hybrid Cloud: This cloud infrastructure is a composite of two or more clouds (private, community, or public) that remain their own entities but are bound together by standard or proprietary technologies that allow portability of data and applications (e.g., cloud bursting for load balancing between clouds).

[0055] Cloud computing environments are service-oriented with an emphasis on statelessness, low coupling, modularity, and semantic interoperability. At the core of cloud computing is an infrastructure that includes a network of interconnected nodes.

[0056] 1B, an exemplary cloud computing environment 52 is shown. As shown, the cloud computing environment 52 comprises one or more cloud computing nodes 50 with which local computing devices used by cloud consumers, such as, for example, a personal digital assistant (PDA) or mobile phone 54A, a desktop computer 54B, a laptop computer 54C, and / or an automobile computer system 54N, may communicate. The nodes 50 may communicate with each other. The nodes may be physically or virtually grouped in one or more networks (not shown), such as a private cloud, a community cloud, a public cloud, or a hybrid cloud, or combinations thereof, as described herein above. This enables the cloud computing environment 52 to provide infrastructure, platform, and / or software as a service for which the cloud consumer does not need to maintain resources on a local computing device. It will be understood that the types of computing devices 54A-N shown in FIG. 1B are intended as examples only, and that the computing nodes 50 and cloud computing environment 52 can communicate with any type of computerized device via any type of network and / or network-addressable connection (e.g., using a web browser).

[0057] Referring now to Figure 1C, a set of functional abstraction layers provided by cloud computing environment 52 (Figure 1B) is shown. It should be understood in advance that the components, layers, and functions shown in Figure 1C are intended to be illustrative only, and embodiments of the present invention are not limited thereto. As shown, the following layers and corresponding functions are provided:

[0058] Hardware and software layer 60 comprises hardware and software components. Examples of hardware components include: mainframe 61; RISC (reduced instruction set computer) architecture based servers 62; servers 63; blade servers 64; storage devices 65; and networks and networking components 66. In some embodiments, software components include network application server software 67 and database software 68.

[0059] The virtualization layer 70 provides an abstraction layer over which the following examples of virtual entities can be provided: virtual servers 71; virtual storage 72; virtual networks, including virtual private networks 73; virtual applications and operating systems 74; and virtual clients 75.

[0060] In one example, management layer 80 may provide the functions described below. Resource provisioning 81 provides dynamic procurement of computing and other resources utilized to execute tasks within the cloud computing environment. Metering and pricing 82 provides cost tracking as resources are utilized within the cloud computing environment and charging or billing for the consumption of these resources. In one example, these resources may include application software licenses. Security provides identity verification for cloud consumers and tasks, and protection for data and other resources. User portal 83 provides access to the cloud computing environment for consumers and system administrators. Service level management 84 provides cloud computing resource allocation and management so that required service levels are met. Service level agreement (SLA) planning and fulfillment 85 provides advance arrangement and procurement of cloud computing resources where future requirements are anticipated according to SLAs.

[0061] Workload tier 90 provides examples of functions for which a cloud computing environment can be utilized. Examples of workloads and functions that can be provided from this tier include: mapping and navigation 91; software development and lifecycle management 92; virtual classroom instructional delivery 93; data analytics processing 94; transaction processing 95; and mobile desktop 96.

[0062] Any of the nodes 50 in the computing environment 52 and computing devices 54A-N may be a DPS 10.

[0063] As discussed in more detail herein, it is contemplated that some or all of the operations of some of the method embodiments described herein may be performed in an alternative order or not at all, and further, multiple operations may be performed simultaneously or as an internal part of a larger process. Determination and use of eigenvalues ​​and eigenvectors of large scale systems.

[0064] Described herein is an engineering analysis system that generates solutions that can be subsequently applied to solve real-world engineering problems. Such an engineering analysis system may be, for example, a spectral clustering analysis system, a frequency response analysis system, or any other system that utilizes matrices as described herein. More generally, the system described herein can be viewed as a computationally efficient technique for generating spectral embeddings. These embeddings can be used for dimensionality reduction purposes. Spectral Clustering Analysis

[0065] Spectral clustering is a technique that uses roots in graph theory, where it is used to identify communities of nodes in a graph based on the edges that connect them. Clustering is one of the main tasks in unsupervised machine learning. Its goal is to assign unlabeled data to groups with the aim that similar data points are assigned to the same group. Spectral clustering is a popular variant of clustering that uses information from the eigenvalues ​​(spectrum) of a special matrix constructed from a graph or dataset. In particular, spectral clustering applies a clustering algorithm to identify the different communities in which the data resides, but does so on an embedded or sketched version of the initial data. This step, which can also be seen as an idea of ​​dimensionality reduction, is accomplished through the computation of a small number of eigenvectors of the similarity matrix.

[0066] In multivariate statistics and clustering of data, the spectral clustering technique uses the spectrum (eigenvalues) of the similarity matrix of the data to perform dimensionality reduction before clustering in fewer dimensions. The similarity matrix is ​​provided as input and consists of a quantitative assessment of the relative similarity of each pair of points in the data set. In its application to image segmentation, spectral clustering is known as segmentation-based object classification.

[0067] Given an enumerated set of data points, the similarity matrix is ​​a symmetric matrix

number

number

[0068] Spectral clustering may involve partitioning a spring-mass system, where each mass is associated with a data point and each spring stiffness corresponds to an edge weight that describes the similarity between two related data points. Specifically, the eigenvalue problem describing the lateral vibration modes of a spring-mass system is given by:

number

number

[0069] Masses that are tightly connected by springs in a spring-mass system move together from their equilibrium positions 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, which is commonly used for image segmentation. It classifies points into two sets (B 1,B 2 ) into

number

[0070] The mathematically equivalent algorithm is to use the random walk normalized adjacency matrix P=D -1 Take the eigenvector corresponding to the largest eigenvalue of A.

[0071] Knowing the eigenvectors, partitioning is done by computing the median m of the components of the second smallest eigenvector v, and then partitioning all points in v whose components are greater than m into B. 1 and the rest in B 2 This can be done in a variety of ways, such as by placing the subsets in a partition of . The algorithm can be used for hierarchical clustering by repeatedly partitioning the subsets in this manner.

[0072] Spectral clustering is computationally expensive unless fast techniques for computing a large number of eigenvectors (i.e., spectral embedding) are available. State-of-the-art methods rely on the application of the (restart-free) Lanczos method. Their cost scales at best linearly with the input size, but with a doubling factor, which depends quadratically on the number of eigenvectors sought. This becomes a major bottleneck when spectral clustering has to be applied to large datasets. Frequency Response Analysis

[0073] Frequency response analysis is the quantitative measurement 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 the magnitude and phase of the output as a function of frequency compared to the input. Simply put, if a sine wave is injected into a system at a given frequency, a linear system will respond at that same frequency with a specific magnitude and a specific phase angle to the input. Sine wave test signals can be used to measure points on the frequency response of a transfer function or impedance function.

[0074] Estimating the frequency response for a physical system generally involves exciting the system with an input signal, measuring both the input and output time histories, and comparing the two through a process such as a Fast Fourier Transform (FFT). The frequency content of the input signal must encompass the frequency range of interest, because the results are not valid for parts of the frequency range that are not encompassed.

[0075] The frequency response of a system can be measured by applying test signals, for example: · Applying an impulse to a system and measuring its response (looking at the 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 a multi-frequency signal (Non-orthogonal frequency discrete multiplexing of signals (N-OFDM or the same as SEFDM) and OFDM) digitally generated maximal length sequence noise, or its equivalent analog filtered white noise such as pink noise) and calculating the impulse response by deconvolution of this input and output signal of the system.

[0076] Frequency response is typically characterized by the magnitude of the system's response, measured in decibels (dB) or as decimals, and the phase, measured in radians or degrees, versus frequency in radians / second or Hertz (Hz). These response measurements can be plotted in three ways: by plotting the magnitude and phase measurements against two rectangular plots as a function of frequency to obtain a Bode plot; by plotting the magnitude and phase angle on a single pole plot with frequency as a parameter to obtain a Nyquist plot; or by plotting the 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 vs. frequency portion of the Bode plot may be all that is of interest. For control system design, 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, assuming a phase vs. frequency plot is included for the Bode analysis. The form of the frequency response for a digital system (such as an exemplary FFT filter) is periodic with multiple main lobes and side lobes.

[0077] When the system under investigation is nonlinear, applying a purely linear frequency domain analysis does 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 users to analyze complex nonlinear dynamic effects. Nonlinear frequency response methods reveal complex resonances, intermodulation, and energy transfer effects that cannot be seen using purely linear analysis and are becoming increasingly important in the nonlinear world.

[0078] One important field where frequency response analysis holds a promising role is structural mechanics, where complex structures are modeled through partial differential equations. Due to the high cost of traditional approaches, which stems in part from the large frequency interval, it is essential to use techniques that are more efficient in terms of memory consumption and wall-clock time requirements. One such approach is to instead perform a model analysis, where the entire structure is projected (embedded) into a lower dimensional subspace that captures most of the energy of the problem. For example, an object may vibrate mostly along modes corresponding to a few lowest frequencies, thus justifying the calculation of only the eigenvectors associated with the lowest eigenvalues ​​of the matrix pencil associated with the finite element discretization of the complex structure. Calculating Eigenvectors

[0079] Both spectral clustering and frequency response analysis systems (as well as other analysis systems) rely on determining and applying eigenvectors (outputs) of large matrix pencils (inputs) that routinely appear in these types of analysis systems. Described herein is a system for determining and applying such eigenvalues / eigenvectors. The disclosed system allows for practical computation and application of eigenvalues ​​of matrices that may be too large to fit within system memory.

[0080] 2A is a block diagram illustrating the general concept of such a system 200 as described above. The technical analysis system 210 described herein may include a processor 212 and memory 214 that may be executed on the DPS 10 and used by an eigenvalue / eigenvector determiner 230 to perform the various functions described herein. The technical analysis system 210 may utilize as its input a large sparse symmetric matrix pencil(A,M) 205 generated, for example, in a spectral clustering analysis system or a frequency response analysis system, as described above. Solving the eigenvalues, eigenvectors, or eigenpairs of a matrix system allows for application to technical application systems 250 that use, for example, spectral clustering analysis, as may be used, for example, in AI classification applications, or that use, for example, responses to impulse signals, multi-frequency signals (Non-orthogonal frequency discrete multiplexing of signals (N-OFDM or the same as SEFDM) and OFDM), impulse responses by deconvolution of input and output signals of a system, and frequency response analysis, as may be used in audio analysis systems, and building structural analysis systems.

[0081] The symmetric generalized eigenvalue problem is formally defined as follows: Ax = λMx

[0082] In general, for a matrix A, if there exists a vector x that is not all zero and a scalar λ such that Ax=λMx, then x is said to be an eigenvector of A with the corresponding eigenvalue λ. Matrix A can be thought of as a function that maps vectors to new vectors. Most vectors will be completely different at some point when A is applied to them, but the eigenvectors will only change in magnitude. If one draws a line through the origin and the eigenvectors, after mapping, the eigenvectors will still be tangent to the line. The amount by which the vectors are scaled along the line depends on λ. Eigenvectors are useful in describing the dynamics of the system represented by the matrix, including spectral clustering and frequency response analysis.

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

[0084] Because spectral clustering requires the computation of eigenvectors associated with a small number of smallest eigenvalues ​​of the graph Laplacian matrix, the eigenvector computation task can be quite costly in terms of resources such as processor bandwidth and memory, given the large size of modern data collections.

[0085] In addition, due to the large size of some of the structures subjected to frequency response analysis, solving the equations by conventional techniques is extremely costly. An alternative is the modal superposition method, which requires the calculation of eigenvectors associated with a large number of the smallest eigenvalues.

[0086] The present system may be used to determine eigenvector derivatives of Schur complements in a novel manner, known techniques that, more importantly, give less accurate results and also require more memory than the techniques disclosed herein. Matrix Pencil

[0087] In linear algebra,

number

number

number

[0088] A particular case is the linear matrix pencil

number

[0089] A pencil is called regular if there exists at least one value of λ such that det(A-λB)≠0. The eigenvalues ​​of a matrix pencil(A,B) are all complex values ​​λ for which det(A-λB)=0. The set of eigenvalues ​​is called the spectrum of the pencil and is written as σ(A,B). Moreover, a pencil is said to have one or more eigenvalues ​​at infinity if B has one or more zero eigenvalues.

[0090] This disclosure considers a Rayleigh-Ritz projection scheme to determine a selected number of eigenvalues ​​(and, optionally, associated eigenvectors) of a large sparse symmetric pencil. This represents an improvement over the traditional shift-and-flip Lanczos method previously used. The focus is on determining a small number of eigenvalues ​​located immediately to the right (or left) of some real scalar, and the proposed technique is 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.

[0091] Determining the eigenvalues ​​that are located closest to a given scalar is typically accomplished by enhancing a selected projection method with shifts and inversions. When less accuracy is required in approximating the sought eigenvalues, an alternative to shifts and inversions is the Automated Multilevel Substructuring (AMLS) technique. AMLS builds large projected subspaces while avoiding excessive orthogonalization costs, and can be significantly faster than shift-and-invert Krylov subspace techniques in applications where a large number of eigenvalues ​​are sought. Process Flow

[0092] FIG. 3 is a flow chart illustrating a process 300 for finding an eigenpair (eigenvalue, eigenvector) according to some embodiments.

[0093] In operation 305, the engineering analysis system 210 may load / assemble a matrix system, for example, a large sparse symmetric matrix pencil(A,M) 205, which constitutes a domain. The focus here is on assuming that matrices A and M are symmetric and that M is a symmetric positive definite (SPD) matrix. The matrix pencil(A,M) 205 is now denoted by (λ i ,x (i) ), with n eigenpairs denoted by i=1,...,n, where λ i constitute the eigenvalues, and x (i) constitute the eigenvectors. We are interested in the n ev The determination of eigenpairs, which can be expressed as some scalar

number

number

[0094] The loading / assembly of the matrix system is done by using matrices A and M, a scalar n ev , and α, a scalar value from which eigenvalues ​​are computed. It also includes inputting a number of partitions p, and a number of local modes to compute from each substructure (partition).

[0095] The standard technique for computing a small number of eigenpairs of a sparse symmetric matrix pencil is

number

number

number

number

[0096] matrix

number

number

number

[0097] Ax using block format (i) = λ i Mx (i) gives

number

number

number

[0098] From the domain segmentation point of view, the ideal choice is

number

number

[0099] The main goal of an algebraic domain decomposition eigensolver is to solve a subspace

number

[0100] 2D, the first operation performed on matrix pencil (A,M) 205 is to have eigenvalue / eigenvector determiner 230 perform partitioning into sub-domains and reordering the matrix in operation 310. The graph associated with a given matrix pencil may be partitioned into multiple non-overlapping sub-domains with the aid of algebraic graph partitioner 232. The projected subspace may have boundary variables 234 and internal variables 236.

[0101] A practical application of the algebraic domain decomposition eigensolver is the problem of solving the matrix B-λM B and E-λM E The method relies on relabeling the unknowns / formulas of the eigenvalue problem Ax=λMx such that |A|+|M| is block diagonal. This can be easily achieved by applying a graph partitioner to the adjacency graph of the matrices |A|+|M|. In this disclosure, only p-way partitionings are considered, and the partitioner divides the graph into p>1 non-overlapping subregions. The rows / columns of the matrices A and M are then reordered such that the unknowns / formulas associated with the internal variables (i.e., the nodes of the adjacency graph that are only connected to nodes located in the same partition) are listed before the unknowns / formulas associated with the boundary variables (i.e., the nodes of the adjacency graph that are connected to nodes located in neighboring partitions). The sought eigenpairs of the original pencil are then approximated by a Rayleigh-Ritz projection onto the subspace formed by the direct sum of the subspaces associated with the two different types of variables. algebraic domain decomposition

[0102] This disclosure focuses on algebraic domain decomposition eigensolvers, where the concept of domain decomposition is applied directly to the eigenvalue equations. The main challenge of algebraic domain decomposition eigensolvers is the construction of the boundary projected subspace, due to the nonlinear nature of the boundary matrix operator, also known as the "spectral Schur complement".

[0103] In linear algebra and matrix theory, the Schur complement of a block matrix is ​​defined as follows: Let A, B, C, and D be p×p, p×q, q×p, and q×q matrices, respectively, and assume that D is invertible. Such that M is a (p+q)×(p+q) matrix,

number

[0104] Next, the Schur complement of block D of matrix M is M / D:=A-BD -1 C If A is invertible, then the Schur complement of a block A of matrix M is M / A:=C-CA -1 B If A or D are singular, substituting the generalized inverse for the inverse for M / A and M / D yields the generalized Schur complement.

[0105] AMLS, an alternative technique to the one focused on here, can be seen as constructing this subspace by computing a small number of eigenvectors of a linear generalized eigenvalue problem derived by first-order truncation of a Taylor series of the spectral Schur complement expanded around some real shift. As a result, the accuracy provided by AMLS can degrade significantly as one moves away from the selected shift. A proposed alternative is to form a boundary projected subspace by considering a complex rationalization transformation to the original pencil that attenuates the boundary components of those eigenvectors associated with the undesired eigenvalues. While these techniques may indeed result in improved accuracy, they also require multiple matrix decompositions computed in complex arithmetic, which is processor intensive.

[0106] The algorithm disclosed herein: a) preserves the advantages of algebraic domain decomposition eigensolvers, such as reduced orthogonalization cost and inherent parallelism, while at the same time, b) may enhance their accuracy without considering more than one shift. Numerical experiments have demonstrated the efficiency of the proposed technique on serial / distributed memory architectures and its competitiveness against schemes such as shift-and-reverse Lanczos and AMLS combined with p-vertex-based partitioning. Method characteristics

[0107] The main properties of the method include: 1) zeroth-order truncation of the boundary matrix operator; 2) exploitation of Taylor series expansion of the boundary eigenvectors; and 3) reduced orthogonalization cost and improved parallelism, which are summarized below. Zeroth-order truncation of boundary matrix operators

[0108] In contrast to AMLS, the algorithm disclosed herein only considers the zeroth order truncation of the spectral Schur complement (the boundary matrix operator), i.e., only a partial solution of the standard eigenvalue problem is required to generate the boundary projected subspace 234. This approach avoids the need to compute / apply derivatives of the spectral Schur complement, as AMLS does. Use of Taylor series expansion of boundary eigenvectors

[0109] The accuracy of the boundary projected subspace is improved by expanding the (analytical) eigenvectors of the spectral Schur complement through their Taylor series and injecting a small number of leading eigenvector derivatives into the subspace. It is theoretically shown and experimentally verified that injecting up to the second leading eigenvector derivatives results in up to fourth order eigenvalue approximation errors. These eigenvector derivatives are computed independently of each other by utilizing a deflationary Krylov subspace solver. Some details are discussed and constraints on the effective condition number are given. Reduced orthogonalization cost and increased parallelism

[0110] As with domain decomposition schemes such as AMLS and RF-DDES, orthonormalization is applied only to vectors whose length is equal to the number of boundary variables instead of the entire set of domain variables (e.g., as with shift-and-reverse Lanczos). This becomes important when a large number of eigenvalues ​​are to be found and the number of boundary variables is much smaller compared to the total number of equations / unknowns.

[0111] In addition, the projected subspaces associated with the internal variables of each subdomain are constructed independently of each other (hence trivially in parallel) by computing local eigenmodes and then computing up to second-order resolvent expansions. Details on the expected accuracy of the eigenvalue approximations are presented.

[0112] In operation 310, the operation creates a projection of the subspace associated with the boundary variables 234. Typically, this is done by locating the bottleneck, i.e.

number

number

number

[0113] This approach also has the advantages of parallelism and orthogonalization cost, is derivative-free, and may result in higher accuracy, but may or may not reduce computational cost.

[0114] In operation 315, a subspace projection associated with the boundary variables α and n ev As an input, operation 315 has two main suboperations: ev The smallest positive eigenvalues ​​μ(α) and associated eigenvectors y(α) (these eigenpairs are (μ ψ (α),y ψ (α), i=1,...,n ev (1) solving S(α)y(α)=μ(α)y(α) for

number

[0115]

number

number

[0116] FIG. 2E is a graph 240 showing elements within suboperation (1) of operation 315 that utilizes the truncated spectral Schur complement. j Let (ζ) denote the j-th eigenvalue of the matrix, where:

number

[0117] Here, μ j Assume that (λ)=0, where j is known and σ∈[λ - ,λ + ] is μ j (.) containing only the λ root

number

number

[0118] Through its Taylor series around σ, j Expanding (σ) gives:

number

number

[0119] Referring to graph 240' in FIG. 2F, the above discussion assumes that: a) λ is known; b) μ j (λ) = 0; and c) y, dy, d 2 We omit some practical details, including requiring a formula that approximates y, ... In practice, we expect that continuous eigenvalue curves will cross the real axis at continuous λ∈Λ(A,M). Therefore, the algorithm is ev The n smallest nonnegative eigenvalues ​​associated with ev Only eigenvalue curves need to be calculated.

[0120] Operation 320 creates a projection of the subspace associated with the internal variables. The main idea here is to

number

number

number

[0121] This operation takes as input σ := α, κ, and Y := the orthonormal basis returned by the previous (boundary variables) algorithm. It then loops for j=1,...,p:

number

number

[0122] Next, in operation 325, a projection matrix Z is constructed for the Rayleigh-Ritz eigenvalue problem as follows: The matrix Y returned from operation 315 is distributed among the p subdomains,

number

number

number

number

[0123] The projection matrix z can then be written as:

number

[0124] The total memory overhead associated with the jth subregion is κd j +(3d j +s j ) η is equal to the total memory overhead of storing a floating-point scalar. T A-Z, Z T The dimension of MZ) is equal to κp+3η and can be solved by suitable routines in LAPACK, e.g., dsygv.

number

number

[0125] Finally, in operation 330, the Rayleigh-Ritz eigenvalue problem is solved for n ev Individual demands

number

number

[0126] The following description compares the algorithm disclosed herein with the shift and inverse Lanczos when it is applied directly to the pencil (A-αM,M).

[0127] The first difference between these two techniques is the orthogonalization cost. Applying k stages of shift and inversion Lanczos to matrix pencils (S(α),I) and (A-αM,M), respectively, costs O(k 2 s) and O(k 2 Therefore, the algorithm reduces the orthogonalization cost by a factor of n / s, and this difference is n ev This becomes more pronounced as k increases (k ≥ n ev The second difference between our algorithm and the shift and inverse Lanczos algorithm is that the matrix B αis the number of linear system solutions using B. It is straightforward to verify that applying k steps of shift-and-flip Lanczos to the pencil (A-αM,M) requires 2k such linear system solutions. In contrast, the algorithm α Using 3η(M E = 0, 2η) linear solutions are required. Nevertheless, the matrix B α These 3η linear solutions using θ can be performed simultaneously since all right-hand sides are simultaneously available, and therefore the associated cost may be slightly higher than the cost of solving a small number of right-hand sides.

[0128] For the above cost, this algorithm uses block diagonal pencils (B α ,M B ) also adds to the computational cost required to compute the eigenvectors of A-αM,M. On the other hand, the accuracy achieved by shift and flip Lanczos applied directly to the pencil (A-αM,M) can be significantly higher than that of the present algorithm. Nevertheless, in many applications the eigenpairs sought do not require very high accuracy at most, e.g. eigenvalue problems resulting from the discretization need to be solved at most up to the associated discretization error.

[0129] The algorithm is based on domain decomposition and is therefore 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.

[0130] Block diagonal matrix

number

number

number

[0131] In contrast to act 320, act 312 involves a computation with the Schur complement matrix S(α), which is distributed among the p subregions. In this case, point-to-point communication between neighboring subregions is required. Finally, the Rayleigh-Ritz eigenvalue problem is typically small enough that it can be replicated and solved redundantly in all MPI processes.

[0132] The detailed application of the eigenvectors of pencil (A,M) computed by the present system on the spectral clustering is as follows: Let A denote the graph Laplacian matrix of an undirected graph, and let M be a diagonal matrix whose ith diagonal entry holds the degree of the ith vertex of the graph, where x (1) , x (2) , ..., x (k) Let denote the eigenvectors associated with the k algebraically smallest eigenvalues ​​of pencil (A,M), and the columns of matrix X denote the k computed eigenvectors x (1) , x (2) , ..., x (k) where x i Let denote the i-th row of the matrix X. The final step of spectral clustering is to invoke a clustering algorithm (e.g., k-means) to partition the k-dimensional points determined by the rows of the matrix X into a user-given number of clusters. Technical Applications

[0133] One or more embodiments disclosed herein thus provide improvements to the computer art, for example, analytical computers capable of efficiently solving eigenvalue / eigenvector problems for large sparse symmetric matrix pencils, providing technical advantages in terms of computer resources such as processor cycles and memory. Computer-Readable Medium

[0134] The present invention may be a system, method, and / or computer readable medium at any possible level of technical detail of integration. A computer program product may include a computer readable storage medium (or media) having computer readable program instructions that cause a processor to perform aspects of the present invention.

[0135] A computer readable storage medium may be a tangible device that can hold and store instructions for use by an instruction execution device. A computer readable storage medium may be, for example, but not limited to, an electronic storage device, a magnetic storage device, an optical storage device, an electromagnetic storage device, a semiconductor storage device, or any suitable combination of the foregoing. A non-exhaustive list of more specific examples of computer readable storage media includes the following: portable computer diskettes, hard disks, random access memories (RAMs), read-only memories (ROMs), erasable programmable read-only memories (EPROMs or flash memories), static random access memories (SRAMs), portable compact disk read-only memories (CD-ROMs), digital versatile disks (DVDs), memory sticks, floppy disks, mechanically encoded devices such as punch cards or ridge structures in grooves that record instructions, and any suitable combination of the foregoing. Computer-readable storage medium, as used herein, should not be construed as a transitory signal per se, such as electric waves or other freely propagating electromagnetic waves, electromagnetic waves propagating through a waveguide or other transmission medium (e.g., light pulses passing through a fiber optic cable), or electrical signals transmitted over wires.

[0136] The computer-readable program instructions described herein can be downloaded from a computer-readable storage medium to the respective computing / processing device or to an external computer or storage device via a network, such as the Internet, a local area network, a wide area network, and / or a wireless network. The network may include copper transmission cables, optical transmission fiber, wireless transmission, routers, firewalls, switches, gateway computers, and / or edge servers. A network adapter card or network interface in each computing / processing device receives the computer-readable program instructions from the network and transfers the computer-readable program instructions to a computer-readable storage medium in the respective computing / processing device for storage.

[0137] The computer readable program instructions for carrying out the operations of the present invention may be either assembler instructions, instruction set architecture (ISA) instructions, machine instructions, machine dependent instructions, microcode, firmware instructions, state setting data, configuration data for an integrated circuit, or source or object code written in any combination of one or more programming languages, including object oriented programming languages ​​such as Smalltalk, C++, etc., and procedural programming languages ​​such as the "C" programming language or similar programming languages. The computer readable program instructions may be executed entirely on the user's computer, partially on the user's computer as a stand-alone software package, partially on the user's computer and partially on a remote computer, or entirely on a remote computer or server. In the latter scenario, the remote computer may be connected to the user's computer via any type of network, including a local area network (LAN) or a wide area network (WAN), or the connection may be made to an external computer (e.g., via the Internet using an Internet Service Provider). In some embodiments, electronic circuitry including, for example, a programmable logic circuit, a field programmable gate array (FPGA), or a programmable logic array (PLA), may execute computer readable program instructions to personalize the electronic circuitry by utilizing state information of the computer readable program instructions to perform aspects of the invention.

[0138] Aspects of the present invention are described herein with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer readable program instructions.

[0139] These computer readable program instructions may be provided to a processor of a general purpose computer, special purpose computer, or other programmable data processing apparatus to produce a machine such that the instructions, executed via the processor of the computer or other programmable data processing apparatus, create means for implementing the functions / acts specified in one or more blocks of the flowcharts and / or block diagrams. These computer readable program instructions may also be stored on a computer readable storage medium such that the instructions can instruct a computer, programmable data processing apparatus, and / or other device to function in a particular manner such that the computer readable storage medium having the instructions stored thereon includes an article of manufacture including instructions that implement aspects of the functions / acts specified in one or more blocks of the flowcharts and / or block diagrams.

[0140] The computer-readable program instructions may also be loaded into a computer, other programmable data processing apparatus, or other device to cause a series of operational steps to be executed on the computer, other programmable apparatus, or other device to generate a computer-implemented process, whereby the instructions executing on the computer, other programmable apparatus, or other device implement the functions / operations specified in one or more blocks of the flowcharts and / or block diagrams.

[0141] The flowcharts and block diagrams in the drawings illustrate the architecture, functionality, and operation of possible implementations of systems, methods, and computer program products according to various embodiments of the present invention. In this regard, each block in the flowcharts or block diagrams may represent a module, segment, or portion of instructions that includes one or more executable instructions that implement the specified logical function. In some alternative implementations, the functions noted in the blocks may occur in an order different from the order noted in the drawings. For example, two blocks shown in succession may in fact be executed substantially simultaneously, or the blocks may be executed in the reverse order depending on the functionality involved. It should also be noted that each block of the block diagrams and / or flowchart diagrams, as well as combinations of blocks in the block diagrams and / or flowchart diagrams, may be implemented by a dedicated hardware-based system that executes the specified functions or operations, or executes a combination of dedicated hardware and computer instructions. 。 [Item 1] 1. A method for solving eigenvalues ​​of a matrix system utilizing a computing device, comprising: Using a processor of the computing device: loading the matrix system into a memory of the computing device; performing a partitioning into non-overlapping sub-regions and a rearrangement of the matrices of the matrix system; creating a boundary variable projection subspace associated with boundary variables of an adjacency graph of the matrix system, where the boundary variables are related to nodes of the adjacency graph that are connected to nodes located in neighboring partitions; creating an internal variable projection subspace associated with internal variables of the adjacency graph of the matrix system, where the internal variables are related to nodes of the adjacency graph that are connected only to nodes located in the same partition; constructing a projection matrix based on the boundary variable projected subspace and the internal variable projected subspace; determining eigenvalues ​​that solve a Rayleigh-Ritz eigenvalue problem using the projection matrix; and Utilizing said determined eigenvalues ​​in a technical application system. A method comprising: [Item 2] 2. The method according to claim 1, wherein the technical application system is selected from the group consisting of a clustering analysis system and a frequency response analysis system. [Item 3] 3. The method of claim 2, wherein the clustering analysis system is an artificial intelligence (AI) classification system and the frequency response analysis system is selected from the group consisting of an audio analysis system and a building structural analysis system. [Item 4] 2. The method of claim 1, wherein creating the boundary variable projected subspace associated with the boundary variables utilizes a zeroth order truncation of the spectral Schur complement matrix and avoids computing derivatives of the spectral Schur complement matrix. [Item 5] Creating the boundary variable projected subspace associated with the boundary variables comprises: expanding the analytical eigenvectors of the spectral Schur complement matrix through their Taylor series; and injecting the leading eigenvector derivatives into the boundary projected subspace. 2. The method according to item 1, comprising the steps of: [Item 6] 2. The method of claim 1, wherein the matrix system is a large sparse symmetric matrix pencil (A,M) that includes a matrix A and a matrix M, the matrix pencil constituting a domain of the matrix system. [Item 7] The matrices A and M are symmetric, The matrix M is a symmetric positive definite matrix. The method according to item 6. [Item 8] 7. The method of claim 6, wherein the step of performing the region partitioning utilizes an algebraic graph partitioner on an adjacency graph associated with the matrix system. [Item 9] determining an eigenvector associated with each eigenvalue to form an eigenpair; and Utilizing said eigenpairs in said technical application system. 2. The method of claim 1, further comprising: [Item 10] An apparatus for an engineering analysis system comprising: memory; and A processor comprising: loading a matrix system into said memory; performing a partition into non-overlapping sub-regions and a matrix reordering of the matrix system; creating a boundary variable projection subspace associated with boundary variables of an adjacency graph of the matrix system, where the boundary variables are related to nodes of the adjacency graph that are connected to nodes located in neighboring partitions; creating an internal variable projection subspace associated with internal variables of the adjacency graph of the matrix system, where the internal variables are related to nodes of the adjacency graph that are connected only to nodes located in the same partition; constructing a projection matrix based on the boundary variable projected subspace and the internal variable projected subspace; determining eigenvalues ​​that solve a Rayleigh-Ritz eigenvalue problem using the projection matrix; and Utilizing said determined eigenvalues ​​in a technical application system. An apparatus comprising a processor configured to: [Item 11] Item 11. The apparatus of item 10, wherein the technical application system is selected from the group consisting of a clustering analysis system and a frequency response analysis system. [Item 12] Item 12. The apparatus of item 11, wherein the clustering analysis system is an artificial intelligence (AI) classification system and the frequency response analysis system is selected from the group consisting of an audio analysis system and an architectural structural analysis system. [Item 13] Item 11. The apparatus of item 10, wherein creating the boundary variable projected subspace associated with the boundary variables utilizes a zeroth order truncation of the spectral Schur complement matrix and avoids computing derivatives of the spectral Schur complement matrix. [Item 14] Creating a boundary variable projected subspace associated with the boundary variables includes: expanding the analytical eigenvectors of said spectral Schur complement matrix through their Taylor series; and injecting the leading eigenvector derivatives into said boundary projected subspace; Item 11. The device according to item 10, comprising: [Item 15] Item 11. The apparatus of item 10, wherein the matrix system is a large sparse symmetric matrix pencil (A,M) that includes a matrix A and a matrix M, the matrix pencil constituting a domain of the matrix system. [Item 16] The matrices A and M are symmetric, The matrix M is a symmetric positive definite matrix. Item 16. The apparatus according to item 15. [Item 17] Item 16. The apparatus of item 15, wherein performing the region partitioning utilizes an algebraic graph partitioner on an adjacency graph associated with the matrix system. [Item 18] The processor, determining an eigenvector associated with each eigenvalue to form an eigenpair; and Utilizing said eigenpairs in said technical application system. Item 11. The apparatus of item 10, further configured to: [Item 19] A computer program product for a technical analysis system comprising: One or more computer readable storage media and program instructions collectively stored on the one or more computer readable storage media, the program instructions comprising: program instructions for loading a matrix system into said memory; program instructions for performing a partition into non-overlapping sub-regions and a rearrangement of the matrices of the matrix system; program instructions for creating a boundary variable projection subspace associated with boundary variables of an adjacency graph of said matrix system, where said boundary variables are related to nodes of said adjacency graph that are connected to nodes located in neighboring partitions; program instructions for creating an internal variable projection subspace associated with internal variables of the adjacency graph of the matrix system, where the internal variables are related to nodes of the adjacency graph that are connected only to nodes located in the same partition; program instructions for constructing a projection matrix based on the boundary variable projected subspace and the internal variable projected subspace; program instructions for determining eigenvalues ​​that solve a Rayleigh-Ritz eigenvalue problem using the projection matrix; and program instructions for utilizing said determined eigenvalues ​​in a technical application system 1. A computer program product comprising: [Item 20] the matrix system is a large sparse symmetric matrix pencil (A,M) that includes a matrix A and a matrix M, the matrix pencil constituting a domain of the matrix system; Matrices A and M are symmetric; The matrix M is a symmetric positive definite matrix; and performing the region partitioning using an algebraic graph partitioner on an adjacency graph associated with the matrix system; 20. The computer program product of item 19.

[0142] Addendum to the Specification

[0143] A domain decomposition Rayleigh-Ritz algorithm for symmetric generalized eigenvalue problems (This version was last compiled on September 24, 2020.)

[0144] VASSILIS KALANTZIS (IBM Research, Thomas J. Watson Research Center, Yorktown Heights, NY 10598, USA. Email:vcal@ibm.com)

[0145] Abstract. This paper proposes an algebraic domain decomposition scheme to compute a selected number of eigenvalues ​​(and, optionally, the associated eigenvectors) of large sparse symmetric matrix pencils. The projected subspaces associated with the boundary variables are constructed by computing a small number of eigenvectors and associated first and second derivatives of a zeroth-order approximation of a nonlinear matrix-valued boundary operator. The projected subspaces associated with the internal variables are constructed independently in each subdomain by exploiting local eigenmodes and matrix resolvent approximations. The sought eigenpairs are then approximated via Rayleigh-Ritz projection. Some theoretical and practical details are discussed, and upper bounds on the approximation error are provided. Numerical experiments demonstrate the efficiency of the proposed technique on serial and distributed memory architectures, and its competitiveness against schemes such as shift-and-reverse Lanczos, and automated multilevel substructuring combined with p-vertex-based partitioning.

[0146] Keywords: symmetric generalized eigenvalue problem, spectral Schur complement, shift and inverse Lanczos, domain decomposition, algebraic substructuring, parallel computing

[0147] AMS subject classification.65F15, 15A18, 65F50

[0148] 1. Introduction. This paper proposes a Rayleigh-Ritz projection scheme based on algebraic domain decomposition to compute the eigenvalues ​​(and, optionally, the associated eigenvectors) of large sparse symmetric matrix pencils. In particular, our focus is on computing a large number of eigenvalues ​​that lie immediately to the right of some real scalar. Eigenvalue problems of this form appear, among others, in applications such as low frequency response analysis [16, 31] and spectral clustering

[39] . The extension to the case where the desired eigenvalues ​​lie immediately to the left of some real scalar is straightforward.

[0149] Computing the eigenvalues ​​located closest to a given scalar is typically achieved by improving the selected projection method by shifting and inverting [13, 41]. When less accuracy is required in approximating the sought eigenvalues, an alternative to the Krylov subspace technique is the automated multilevel substructuring technique (AMLS) [6, 7, 10, 28]. From an algebraic point of view, AMLS performs a Rayleigh-Ritz projection onto a large projected subspace while avoiding excessive orthogonalization costs, and has been demonstrated as a good alternative to the shift-and-invert Krylov subspace technique when a very large number of eigenvalues ​​are sought

[29] . An analysis of the AMLS algorithm for elliptic PDE problems from a transformation point of view can be found in [7], where tools from the mathematical theory of substructured domain decomposition for elliptic PDEs are considered to understand how AMLS scales as the mesh is refined.

[0150] This paper focuses on algebraic domain decomposition eigensolvers, where the concept of domain decomposition is applied directly to the eigenvalue equations. The algebraic domain decomposition eigensolver starts by invoking an algebraic graph partitioner to partition the graph associated with a given matrix pencil into multiple non-overlapping subdomains. The variables in each subdomain are then classified into two different categories: a) variables that are connected to variables located only in the same subdomain (interior), and b) variables that are connected to both local variables and variables located in neighboring subdomains (boundary). The eigenpairs sought are then approximated by Rayleigh-Ritz projections onto the subspace formed by the union of the subspaces associated with the two different types of variables. A thorough analysis of the algebraic domain decomposition eigensolver can be found in

[20] .

[0151] The main challenge of algebraic domain decomposition eigensolvers is the construction of the projected subspace associated with the boundary variables due to the nonlinear nature of the boundary matrix operator, also known as the "spectral Schur complement" [21, 23]. From a purely algebraic point of view, AMLS constructs this subspace by computing a small number of eigenvectors of a linear generalized eigenproblem involving the Schur complement and its first derivatives [6]. As a result, the accuracy provided by AMLS can degrade significantly as we turn our attention to computing eigenvalues ​​located away from the origin. A recently proposed alternative is to form the projected subspace associated with the boundary variables by applying a complex rationalization transformation to the original pencil, so as to vanish the components associated with the undesired eigenvalues. See, for example, the RF-DDES algorithm

[24] and the discussion in [22, 27]. While these techniques can indeed result in improved accuracy, they require multiple matrix decompositions computed in complex arithmetic.

[0152] In this paper, we propose an algorithm that preserves the advantages of algebraic domain decomposition eigensolvers, such as reduced orthogonalization cost and inherent parallelism, while at the same time increasing their accuracy without considering shifts more than one.

[0153] The important properties of the proposed scheme are summarized below.

[0154] 1) Zeroth-order truncation of the boundary matrix operator. In contrast to AMLS, the algorithm proposed in this paper generates the projected subspace associated with the boundary variables by partially solving the standard eigenvalue problem with the Schur complement. This approach avoids the need to compute / apply derivatives of the spectral Schur complement.

[0155] 2) Exploiting Taylor series expansion of boundary eigenvectors. The accuracy of the projected subspace associated with the boundary variables is improved by expanding the (analytical) eigenvectors of the spectral Schur complement through their Taylor series and injecting a small number of leading eigenvector derivatives into the subspace. We show theoretically (and experimentally verify) that injecting at most second order eigenvector derivatives leads to eigenvalue approximations with absolute error upper bounds reduced to fourth order. These eigenvector derivatives are computed independently of each other by exploiting a deflationary Krylov subspace solver.

[0156] 3) Reduced orthogonalization cost and improved parallelism. Similar to domain decomposition schemes such as AMLS and RF-DDES, orthonormalization is applied only to vectors whose length is equal to the number of boundary variables instead of the entire set of domain variables. This becomes especially important when a large number of eigenvalues ​​are sought and the number of boundary variables is much smaller compared to the total number of equations / unknowns. In addition, the projected subspaces associated with the internal variables of each subdomain are constructed independently from each other by computing local eigenmodes and at most quadratic resolvent expansions. We report experiments performed on serial / distributed memory architectures and demonstrate the suitability of the proposed technique in high-performance computing settings.

[0157]

number

[0158] 1.1 Notation and Overview. The eigenvalue problems considered in this paper are of the form Ax = λMx, where the matrix

number

number

[0159] The outline of this paper is as follows: Section 2 gives a brief introduction to Rayleigh-Ritz projected subspaces from the domain decomposition point of view. Section 3 presents a spectral Schur complement technique that approximates one or more eigenpairs of pencil (A,M) independently from each other. Section 4 extends these results to Rayleigh-Ritz projection that approximates all sought eigenpairs from a single projected subspace. Section 5 presents numerical experiments performed on serial and distributed memory machines. Finally, conclusions are given in Section 6.

[0160] 2. Rayleigh-Ritz projection from the domain segmentation perspective. The standard technique for computing a small number of eigenpairs of a sparse symmetric matrix pencil is

number

number

[33] by applying the Rayleigh-Ritz technique to ev The desired eigenvalues

number

number

number

[0161] The matrices A and M are assumed to be partitioned into a 2x2 block format.

number

[0162] Here, B and M B is a square matrix of size d × d, and E and M E is a rectangular matrix of size d×s, and C and M C is a square matrix of size s × s, and n = d + s. Without loss of generality, n ev Assume that s ≦ s. Similarly, the eigenvalue λ i The eigenvector x associated with (i) can be written as follows:

number

[0163] Using the block format above, Ax (i) = λ i Mx (i) Rewriting this gives us:

number

[0164] On the other hand, from the second line of (2.1), u (i) Eliminating

[0165]

number

[0166] From these, λ i and y (i) It is possible to determine x (i) , u (i) The missing d×1 part of is then restored by solving the linear system

[0167]

number

[0168] From the domain segmentation point of view, the ideal choice is

number

number

[0169] The main goal of an algebraic domain decomposition eigensolver is to construct a projected subspace that ideally contains the subspace in (2.3).

[0170] 2.1. Region partition sorting. A practical application of the algebraic domain decomposition eigensolver is the problem of solving the matrix B-λM BWe rely on relabeling the unknowns / equations in the eigenvalue problem Ax = λMx such that is block diagonal. This can be easily achieved by applying a graph partitioner to the adjacency graph of the matrix |A| + |M|, e.g. [25, 34].

[0171] In this paper, we consider only p-partitioning, where the partitioner divides the graph into p>1 non-overlapping subregions. The rows / columns of the matrices A and M are then reordered such that the unknowns / expressions associated with the internal variables (i.e., nodes of the adjacency graph that are only connected to nodes located in the same partition) are listed before the unknowns / expressions associated with the boundary variables (i.e., nodes of the adjacency graph that are connected to nodes located in neighboring partitions). The permuted matrices A and M can be written as follows (the eigenvalues ​​of pencil(A,M) in (2.4) are identical to the eigenvalues ​​before the symmetric permutation. If the eigenvectors are also of interest, then these can be written as P T Throughout the remainder of the paper, we work with the permuted matrix shown in (2.4).

number

[0172] The internal and boundary variables in the jth subregion are denoted by d j and s j In this case, the matrix B j and

number

number

number

number

number

number

number

number

[0173] Using the Schur complement framework, we can write the vector u (i) and y (i) teeth,

number

number

number

number

[0174] Partial Vector

number

number

number

[0175] 3. Scheme for approximating a single eigenpair. This section

number

number

number

number

[0176] 3.1. Spectral Schur Complement We consider a univariate, nonlinear matrix-valued function

number

number

[0177] Definition 3.1.For j=1,...,s, define the following scalar-vector pair:

number

number

number

number

[0178] Each function μ j (ζ) (hereafter referred to as the "eigenvalue curve") is the

number

[0179] Note 1. Throughout this paper, μ j Let μ be the eigenvalue curve whose (0) is equal to the j-th algebraically smallest eigenvalue of the matrix S(0). j We denote this by (ζ). This idea can be extended to real scalars other than the origin.

[0180] 3.2. Behavior of eigenvalue curves at their own poles. In general, (B,M B ) eigenvalues ​​(if any) of μ j It is not possible to determine which poles of (ζ) are in (B,M B ) can be used to determine the number of eigenvalue curves that remain analytical as one approaches the eigenvalues ​​of

[0181] Definition 3.2. Matrix Pencil(B,M B ) is an eigenpair of (δ l ,v (l) ), denoted by l=1,...,d, where V=[v (1) ,v (2) ,...,v (d) ] is V T M B It is scaled so that V=I.

[0182]

number

number

number

number

[0183] ζ→δ k Then, all the eigenvalue curves except the above ones can be analytically expressed as δ k The eigenvalue curves are monotonically decreasing over their domain. Further details on the above discussion can be found in [23, Theorem 4.1] and [30, Theorem 2.1] for M = I. The extension to the case M ≠ I is straightforward.

[0184] 3.3.y j Taylor series expansion of (λ). Following (2.2), the scalar

number

[30] . An illustration of this concept is shown in Figure 3.1. In that case, the corresponding eigenvector y j (λ) is parallel to the lowest s × 1 subvector of eigenvector x.

[0185] Here, the eigenpair (μ j (σ), y j (σ)), where σ∈[λ - ,λ + ], and μ jThe latter interval lies between two consecutive poles of (ζ), if such poles exist. (In fact, the interval [λ - ,λ + μ located inside j The poles of (ζ) pose no threat as long as σ is chosen carefully. Ideally, σ ≡ λ, and therefore y j (σ)=y j (λ)=y. In reality,

number

[0186]

number

number

[0187]

number

[0188] In that case, the eigenvector x can be written as:

number

[0189] The above formula for eigenvector x is y j (σ) as its zeroth approximation, and y evaluated at σ jBy considering higher derivatives of (ζ), y j We propose that an improved approximation of (λ) can be achieved.

[0190] 3.4. Rayleigh quotient approximation of λ

[0191] Definition 3.3.

number

number

[0192] Definition 3.4. SPD n×n matrix M and non-zero vector

number

number

[0193] 3.4.1. Basic approximations.

[0194] Proposition 3.5. λ∈Λ(A,M) is μ j (λ)=0, and λ - ≦σ≦λ≦λ + where λ is in the interval [λ - ,λ + μ located in j is the unique root of (ζ). In addition,

number

number

number

[0195] next,

number

number

[0196] Proof.

number

number

number

[0197] The slope of the Rayleigh quotient is

number

number

[0198] Here, E λ =E σ -(λ-σ)M E Describe the vector

number

[0199] In that case, the difference

number

number

[0200] Where:

number

number

number

number

number

[0201] identity

number

number

number

number

number

[0202] Proposition 3.5 states that in the interval [λ - ,λ + ] is (including multiple eigenvalues) μ j It remains valid even if it contains more than one root of (ζ). However, it can only approximate one of these eigenvalues. Moreover, in the interval [λ - ,λ + ] is μ j It is possible to include poles of (ζ), but σ must be algebraically smaller than these poles.

[0203] 3.4.2.

number

number

number

[0204] More specifically, the approximate eigenvector

number

number

number

number

[0205] Following the reasoning in Proposition 3.5 (omitting the intermediate steps), it can be shown that

number

[0206] Eigenvector approximation shown in (3.5)

number

number

number

[0207] 3.4.3. Reduction of the asymptotic order of the upper bound. The analysis in Proposition 3.5 shows that, regardless of the value of τ, in the limit σ → λ, the term

number

number

number

[0208]

number

number

number

number

number

number

number

number

[0209] Proposition 3.6. λ∈Λ(A,M) is μ j (λ)=0, and λ - ≦σ≦λ≦λ + where λ is in the interval [λ - ,λ + μ located in j is the unique root of (ζ). In addition,

number

number

number

number

[0210] next,

number

number

[0211] Proof. First,

number

[0212] Next, the difference

number

number

[0213]

number

[0214] Following the same reasoning in Proposition 3.5 (omitting the intermediate steps), we get:

number

[0215] Proposition 3.6 states that the resolvent

number

number

[0216] 4. Computation of Multiple Eigenpairs. The technique discussed in Section 3 can be extended to the simultaneous approximation of all desired eigenpairs. More specifically, we use the n eigenpairs immediately to the right of the real scalar α. ev The desired eigenvalues ​​are

number

number

number

number

number

number

number

number

[0217] For the time being, we assume that the eigenvalue curves do not cross each other and set σ = α. Since the eigenvalue curves are monotonically decreasing, the eigenvalue

number

number

number

number

number

[0218]

number

[0219] In reality, the eigenvalue curves that intersect with each other are n of the matrix S(α). ev Each one of the algebraically smallest nonnegative eigenvalues ​​has value

number

[0220] 4.1. The Rayleigh-Ritz Algorithm. The theoretical results presented in Section 3 are used to approximate eigenvectors

number

number

[0221]

number

[0222] Algorithm 4.1 uses a user-supplied scalar

number

number

[0223] 4.2. Construction of projected subspaces associated with boundary variables. Algorithm 4.2 is the construction of n ev We start by computing the eigenvectors associated with the smallest non-negative eigenvalues. These can be computed by any suitable sparse eigensolver, for example the (block) Lanczos method combined with shift-and-inversion

[13] or the (generalized) Davidson method

[36] . The algorithm is based on the derivative

number

[0224]

number

[0225] 4.3. Constructing the projected subspace associated with the internal variables. Algorithm 4.3 provides a formal description of the procedure followed for constructing the projected subspace associated with the internal variables. The routine "

number

number

[40] .

[0226]

number

[0227] 4.3.1. The Rayleigh-Ritz eigenvalue problem. The matrix Y returned by Algorithm 4.2 is distributed among the p subdomains,

number

number

number

number

[0228] The projection matrix Z can then be written as:

number

[0229] The total memory overhead associated with the jth subregion is κd j +(3d j +s j ) η is equal to the total memory overhead of storing a floating-point scalar. T A-Z, Z T The dimension of MZ is equal to κp+3η and can be solved by suitable routines in LAPACK, e.g., dsygv[3].

number

number

[0230] 4.4. Comparison of Algorithm 4.1 with Shift-and-Inverse Lanczos A natural question is how Algorithm 4.1 compares to the shift-and-inverse Lanczos when it is applied directly to the pencil (A-αM,M).

[0231] The first important difference between these two techniques is the orthogonalization cost. Applying k stages of shift and inversion Lanczos to the matrix pencils (S(α),I) and (A-αM,M), respectively, costs O(k 2 s) and O(k 2Therefore, Algorithm 4.1 reduces the orthogonalization cost by a factor of n / s, and this difference is n ev This becomes more pronounced as k increases (k ≥ n ev because.

[0232] The second important difference between Algorithm 4.1 and the shift-and-reverse Lanczos is that we use matrix B α is the number of linear system solutions using B. It is straightforward to verify that applying k steps of shift-and-flip Lanczos to pencil (A-αM,M) requires 2k such linear system solutions. In contrast, Algorithm 4.1 α Using 3η(M E = 0, then 2η) linear solutions are required. Nevertheless, these 3η linear solutions can be executed simultaneously since all right-hand sides are simultaneously available. Therefore, the associated cost can be slightly higher than the cost of solving a small number of right-hand sides. For the above cost, we use block diagonal pencils (B α ,M B ) one must also add the computational cost required to compute the eigenvectors of

[0233] On the other hand, the accuracy achieved by shift-and-reverse Lanczos can be significantly higher than that of Algorithm 4.1. Nevertheless, in many applications the eigenpairs sought do not require very high accuracy at most, e.g., eigenvalue problems resulting from discretization need to be solved at most up to the associated discretization error.

[0234] 4.5. Parallel computing. Algorithm 4.1 is based on domain decomposition and is therefore well suited for execution in modern distributed memory environments: for example, each one of the p subdomains can be mapped to a separate MPI process. matrix

number

number

number

[24] for a similar discussion in the context of domain decomposition eigensolvers, and

[15] for a general discussion on parallel computing and domain decomposition.

[0235] In contrast to Algorithm 4.3, Algorithm 4.2 involves a computation with the Schur complement S(α), which is distributed among the p subdomains. In this case, point-to-point communication between neighboring subdomains is required. Finally, the Rayleigh-Ritz eigenvalue problem is typically small enough that it can be replicated and solved redundantly in all MPI processes.

[0236] 5. Numerical Experiments Our set of experiments is carried out in the Matlab® environment (version R2018b) using 64-bit arithmetic on a single core of a computing system equipped with an Intel Haswell E5-2680v3 processor and 32 GB of system memory.

[0237] The eigenvalue of pencil (A,M) is λ 1 ≦λ 2 <...<λ n and throughout the remainder of this section, n ev algebraically smallest eigenvalues

number

[0238] We consider three different choices for setting span(Y) in Algorithm 4.2:

number

[0239] The eigenpairs of matrix S(α) are at most 1.0×10 -6 The residual norm is calculated up to

[0240] 5.1. Model Problem. Our first test case consists of a five-point stencil finite difference discretization of the Dirichlet eigenvalue problem,

number

[0241] Figure 5.1 shows

number

[0242]

number

[0243]

number

[0244] Table 5.1 shows the ev = 50 and for varying p, we list the max / min relative eigenvalue errors. In summary, increasing the number of boundary variables results in lower accuracy unless more eigenvector derivatives are computed.

[0245] Figure 5.2 shows the case where U is set as in Algorithm 4.3,

number

[0246] 5.2. General Pencil. Throughout the remainder of this section, we assume that the matrices U and Y are set as described in Algorithm 4.1. Details for the numerical approximation of the first and second derivatives of each computed eigenvector of the matrix S(α) are given in the Appendix. (The linear system solver we have chosen has a 1.0×10 -3and a maximum of 5 preconditioned iterations (whichever occurs first). The convergence criterion is applied to the unpreconditioned residuals. For preconditioning, we use a matrix S(α) combined with deflation.

[0247] We consider the application of Algorithm 4.1 to a model problem and a set of matrix pencils obtained by the SuiteSparse (https: / / sparse.tamu / edu / ) Matrix Collection [9] and the Elses (http: / / www.elses.jp / matrix / ) matrix library

[19] .

[0248] [Table 1]

[0249] Details can be found in Table 5.2. Matrix pencil "FEmesh" represents a finite element discretization of the Dirichlet eigenvalue problem for the [-1,1]x[-1,1] plane, where the Laplacian is discretized using linear elements with a target maximum mesh edge length of h=0.05. Matrix "FDmesh" represents a 250x250 finite difference discretization of the Dirichlet eigenvalue problem discussed in the previous section. Except for matrix pencils "FEmesh", "Kuu / Muu", and "qa8fk / qa8fm", all other pencils are in standard form.

[0250] [Table 2]

[0251] [Table 3]

[0252] Table 5.3 shows the n returned by Algorithm 4.1 for all pencils listed in Table 5.2. ev= 50 algebraically smallest eigenvalues. Naturally, increasing κ improves the overall accuracy, and this improvement is even greater as we also improve the boundary projected subspace. However, the choices '{y,dy}' and '{y,dy,d 2 Note that the entries associated with {y} are now closer to each other than those shown in section 5.1. This is due to: a) the inexact calculation of the first and second derivatives, and b) the asymptotic order of the error is the same in both cases. Because the resolvent

number

[0253] Finally, we compare the accuracy achieved by Algorithm 4.1 with our variant of the AMLS method (called “p-way AMLS” and

number

number

number

[0254] 5.3. Comparison to shift and reverse Lanczos. This section compares Algorithm 4.1 (transformation "{y, dy, d 2 We present a real-time comparison between the implicit restarted shift-and-reverse Lanczos with full orthogonalization applied directly to the pencils (A,M). We call the latter “IRSIL”. We define the maximum dimension of the Krylov subspace as 2n ev As our test matrix, we choose a 5-point 506 × 296 finite difference discretization of the Dirichlet eigenvalue problem.

[0255] Figure 5.4 shows the n ev We plot a series of real-time plots of the individual steps of both IRSIL and Algorithm 4.1 when utilizing these schemes to approximate n algebraically smallest eigenvalues. In total, IRSIL takes ev = 50, 100, 150 and n ev = 200 algebraically smallest eigenvalues ​​and a maximum residual norm of 1.0×10 -8 Algorithm 4.1, which uses the default selection p=8, required 4.9, 6.7, 8.7, and 10.5 seconds, respectively, to approximate the associated eigenvectors up to n. ev As n increases, IRSIL becomes increasingly slower than Algorithm 4.1 due to its increasing orthogonalization cost. For example, evFor n = 200, IRSIL spent about 10 seconds in orthogonalization, while algorithm 4.1 required only about a quarter of a second. ev The maximum eigenvalue error returned by Algorithm 4.1 for any tested values ​​is O(10 -4 ) Figure 5.4 also shows that p varies and n ev We plot a series of execution timings of Algorithm 4.1 when p = 200 is fixed (right sub-figure). Increasing the value of p results in a larger number of boundary variables (hence increased orthogonalization cost), while solving the linear system with the matrix S(α) also becomes more expensive. On the other hand, larger values ​​of p generally reduce the amount of time spent in Algorithm 4.3. For p = 64, the maximum error returned by Algorithm 4.1 is of order O(10 -3 Note that in

[0256] Table 5.4 lists the total number of iterations required by Lanczos with implicit restarts: a) to compute the required eigenvectors of matrix S(α), and b) to solve the original eigenvalue problem directly. For Algorithm 4.1, we considered settings p=4, 16, and p=32, respectively. In summary, we observe two patterns. First, as p increases, the number of iterations required by Lanczos also increases. In particular, increasing the size of the Schur complement typically results in a denser spectrum for matrix S(α).

[0257]

number

[0258]

number

[0259] [Table 4]

[0260] The main distributed memory implementation of Algorithm 4.1 is built on top of the PETSc library [5]. (Our code is built on the implementation presented in

[24] and compiled using real arithmetic. The source files were compiled with the Intel MPI compiler mpiicpc using the -O3 optimization level.) The linear system solution using distributed matrices A-αM and S(α) is computed by the Multifrontal Massively Parallel Sparse Direct Solver [2], while the block diagonal matrix B αThe linear system solution using was computed by MKL PARDISO [3].) Message passing between different processes was achieved via the Message Passing Interface (MPI), a standard application programming interface for message passing applications

[14] . We consider only one thread per MPI process, and the number of MPI processes is equal to the number of subdomains p. To allow a fair comparison, we implemented our version of IRSIL instead of the one in the PETSc-based state-of-the-art eigenvalue package SLEPc

[17] . Our experiments were run on the Mesabi cluster at the Minnesota Supercomputing Institute (https: / / www.msi.umn.edu / ).

[0261] Figure 5.6 shows the ev We plot the distributed memory real time of different modules of IRSIL and Algorithm 4.1 for p = 200 and p = 2, 4, 8, 16, 32 and p = 64. The corresponding efficiency for some of these values ​​of p are listed in the attached table.

[0262]

number

[0263]

number

[0264] Finally, Figure 5.7 shows the mesh size

number

[0265]

number

[0266] Moreover, ev Larger values ​​of are more challenging in terms of scalability.

[0267] 6. Summary and Future Work. This paper presents a domain decomposition technique for the computation of a small number of eigenpairs of a symmetric generalized eigenproblem. The proposed technique is based on Rayleigh-Ritz projection into a subspace formed by separating the original eigenproblem into two separate subproblems, each one associated with the boundary and internal variables of the domain, respectively. The part of the subspace associated with the boundary variables of the global domain is formed by the eigenvectors and associated derivatives of the zeroth order approximation of the nonlinear boundary matrix-valued operator. On the other hand, the part of the subspace associated with the internal variables is formed in parallel between different subdomains by exploiting the approximation of the local eigenmodes and resolvent expansions.

[0268] Future work includes investigating the impact that the number of bounding variables has on the accuracy of the techniques proposed in this paper. Parallel implementations using additional levels of parallelism (both distributed and shared memory), possibly combined with harmonic Rayleigh-Ritz projections, are in our near-term plans.

[0269] Acknowledgements. The authors would like to thank the anonymous referees whose detailed comments and suggestions significantly improved the quality of this manuscript. The work discussed in this manuscript was supported by the International Business Machines Corporation through its Herman H. Goldstine Postdoctoral Fellowship program. The authors would like to thank the Minnesota Supercomputing Institute (MSI) at the University of Minnesota for providing resources that contributed to the research results reported in this paper. Additionally, the authors are indebted to Daniel Kressner, Anthony Austin, and Yousef Saad for their feedback and encouragement to pursue the research discussed in this manuscript.

[0270] Appendix A. Analytical formulation.

[0271] In this section, the matrix

number

[0272] Proposition A.1. Matrix

number

number

number

number

number

number

number

number

number

[0273] Normalization condition

number

number

number

number

[0274] NOTE 2 If M is equal to the identity matrix, the first and second derivatives of the matrix-valued function S(ζ) evaluated at σ are given by the block diagonal matrix

number

number

[0275] Appendix B. Calculation of Eigenvector Derivatives.

[0276] Each eigenvector y iThe computation of the first and second derivatives of (σ), i = 1, 2, ..., s, requires the application of an iterative solver, e.g., the MINimum RESidual (MINRES) Krylov subspace method [8, 32], to the solution of a single linear system,

number

number

number

[0277] matrix S(σ)-μ i (σ)I's eigenvalues ​​are {μ k (σ)=μ i (σ)} k=1,2,...,s is equal to.

[0278] n of S(σ) ev The calculated eigenvalues ​​are

number

number

number

number

number

number

number

number

[0279] Details on the deflation MINRES can be found in [11, 12].

[0280] If we already have a symmetric factorization of S(σ) at hand, this can be solved as a linear system (recall that MINRES requires an SPD preconditioner; if S(σ) is not SPD, alternative Kyrlov subspace iterative linear system solvers should be considered).

number

number

[0281] matrix

number

number

number

[0282] System B.1. Find the eigenvalues ​​of matrix S(σ) as follows:

number

[0283] In that case, the matrix

number

number

[0284] Similarly, the preconditioned matrix

number

number

[0285]

number

[0286] n of matrix S(σ) ev If the calculated eigenvectors are not deflated, the denominators in system B.1 are, respectively,

number

number

[18] ; see, for example, [35, 37].

[0287] Figure B.1 shows the maximum eigenvalues ​​of 1.0×10 for a 253×148 finite difference discretization of the Dirichlet eigenvalue problem. -8 Up to a tolerance equal to

number

[0288] References

[0289] [1]Intel (r) fortran compiler xe 14.0 for linux.

[0290] [2]PR AMESTOY, IS DUFF, J.KOSTER, AND J.-Y.L'EXCELLENT, A fully asynchronous multi-frontal solver using distributed dynamic scheduling, SIAM Journal on Matrix Analysis and Applications, 23 (2001),pp. 15-41.

[0291] [3]E. ANDERSON, Z. BAI, C. BISCHOF, S. BLACKFORD, J. DONGARRA, J. DU CROZ, A. GREEN-BAUM, S. HAMMARLING, A. MCKENNEY, AND D. SORENSEN, LAPACK Users' guide, vol. 9, SIAM, 1999.

[0292] [4]A. ANDREW AND R. TAN, Computation of derivatives of repeated eigenvalues and the corre-sponding eigenvectors of symmetric matrix pencils, SIAM Journal on Matrix Analysis and Applications, 20 (1998) , pp. 78-100.

[0293] [5]S. BALAY, S. ABHYANKAR, M. F. ADAMS, J. BROWN, P. BRUNE, K. BUSCHELMAN, L. DALCIN, A. DENER, V. EIJKHOUT, W.D. GROPP, D. KARPEYEV, D. KAUSHIK, M.G. KNEPLEY, D.A. MAY, L.C. MCINNES, R.T. MILLS, T. MUNSON, K. RUPP, P. SANAN, B.F. SMITH, S. ZAMPINI, H. ZHANG, AND H. ZHANG, PETSc users manual, Tech. Rep. ANL-95 / 11-Revision 3.11, Argonne National Laboratory, 2019.

[0294] [6]C. BEKAS AND Y. SAAD, Computation of smallest eigenvalues using spectral Schur comple-ments, SIAM J. Sci. Comput., 27 (2006) , pp. 458-481.

[0295] [7]J.K. BENNIGHOF AND R.B. LEHOUCQ, An automated multilevel substructuring method for eigenspace computation in linear elastodynamics, SIAM J. Sci. Comput., 25 (2004) , pp. 2084-2106.

[0296] [8]S. -C.T. CHOI, C.C. PAIGE, AND M.A. SAUNDERS, MINRES-QLP: A Krylov subspace method for indefinite or singular symmetric systems, SIAM Journal on Scientific Computing, 33 (2011) , pp. 1810-1836.

[0297] [9]T.A. DAVIS AND Y. HU, The University of Florida sparse matrix collection, ACM Trans. Math. Softw., 38 (2011) , pp. 1: 1-1: 25.

[0298]

[10] W. GAO, X.S. LI, C. YANG, AND Z. BAI, An implementation and evaluation of the AMLS method for sparse eigenvalue problems, ACM Trans. Math. Softw., 34 (2008) , pp. 20: 1-20: 28.

[0299]

[11] A. GAUL, M.H. GUTKNECHT, J. LIESEN, AND R. NABBEN, A framework for deflated and augmented Krylov subspace methods, SIAM Journal on Matrix Analysis and Applications, 34 (2013) , pp. 495-518.

[0300]

[12] A. GAUL AND N. Preconditioned recycling Krylon subspace methods for self-adjoint problems, Electronic Transactions on Numerical Analysis, 44 (2015) , pp. 522-547.

[0301]

[13] R.G. GRIMES, J.G. LEWIS, AND H.D. SIMON, A shifted block Lanczos algorithm for solv-ing sparse symmetric generalized eigenproblems, SIAM J. Matrix Anal. Appl., 15 (1994) , pp. 228 272.

[0302]

[14] W. GROPP, W.D. GROPP, E. LUSK, A.D.F.E.E. LUSK, AND A. SKJELLUM, Using MPI: portable parallel programming with the message-passing interface, vol. 1, MIT press, 1999.

[0303]

[15] W.D. GROPP, Parallel computing and domain decomposition, in Fifth International Sympo-sium on Domain Decomposition Methods for Partial Differential Equations, Philadelphia, PA, 1992, pp. 349-361.

[0304]

[16] A. HANNUKAINEN, J. MALINEN, AND A. OJALAMMI, Efficient solution of symmetric eigenvalue problems from families of coupled systems, arXiv preprint arXiv: 1806.07235, (2018) .

[0305]

[17] V. HERNANDEZ, J.E. ROMAN, AND V. VIDAL, Slepc: A scalable and flexible toolkit for the solution of eigenvalue problems, ACM Transactions on Mathematical Software (TOMS) , 31 (2005) , pp. 351-362.

[0306]

[18] M.R. HESTENES AND E. STIEFEL, Methods of conjugate gradients for solving linear systems, Journal of research of the National Bureau of Standards, 49 (1952) , pp. 409-436.

[0307]

[19] T. HOSHI, H. IMACHI, A. KUWATA, K. KAKUDA, T. FUJITA, AND H. MATSUI, Numerical aspect of large-scale electronic state calculation for flexible device material, Japan Journal of Industrial and Applied Mathematics, 36 (2019) , pp. 685-698.

[0308]

[20] V. KALANTZIS, Domain Decomposition Algorithms for the Solution of Sparse Symmetric Gen-eralized Eigenvalue Problems, PhD thesis, University of Minnesota, 2018.

[0309]

[21] V. KALANTZIS, A spectral Newton-Schur algorithm for the solution of symmetric generalized eigenvalue problems, Electron. Trans. Numer. Anal., 52 (2020) , pp. 132-153.

[0310]

[22] V. KALANTZlS, J. KESTYN, E. POLIZZI, AND Y. SAAD, Domain decomposition approaches for accelerating contour integration eigenvalue solvers for symmetric eigenvalue problems, Nu-merical Linear Algebra with Applications, 25 (2018) , p. e2154.

[0311]

[23] V. KALANTZIS, R. LI, AND Y. SAAD, Spectral Schur complement techniques for symmetric eigenvalue problems, Electronic Transactions on Numerical Analysis, 45 (2016) , pp. 305-329.

[0312]

[24] V. KALANTZIS, Y. XI, AND Y. SAAD, Beyond Automated Multilevel Substructuring: domain decomposition with rational filtering, SIAM Journal on Scientific Computing, 40 (2018) , pp. C477-C502.

[0313]

[25] G. KARYPIS AND V. KUMAR, A fast and high quality multilevel scheme for partitioning irregular graphs, SIAM Journal on Scientific Computing, 20 (1998) , pp. 359-392.

[0314]

[26] T. KATO, Perturbation theory for linear operators, 2nd ed., Grundlehren Math. Wiss., Springer, Berlin, 1976.

[0315]

[27] J. KESTYN, V. KALANTZIS, E. POLIZZI, AND Y. SAAD, Pfeast: a high performance sparse eigen-value solver using distributed-memory linear solvers, in SC'16: Proceedings of the Inter-national Conference for High Performance Computing, Networking, Storage and Analysis, IEEE, 2016, pp. 178-189.

[0316]

[28] J.H. Ko AND Z. BAI, High-frequency response analysis via algebraic substructuring, Interna-tional Journal for Numerical Methods in Engineering, 76 (2008) , pp. 295-313.

[0317]

[29] L. KOMZSIK AND T. ROSE, Parallel methods on large-scale structural analysis and physics appli-cations substructuring in MSC / NASTRAN for large scale parallel applications, Computing Systems in Engineering, 2 (1991) , pp. 167 -173.

[0318]

[30] S. LUI, Kron's method for symmetric eigenvalue problems, Journal of computational and applied mathematics, 98 (1998) , pp. 35-48.

[0319]

[31] K. MEERBERGEN AND Z. BAI, The Lanczos method for parameterized symmetric linear systems with multiple right-hand sides, SIAM Journal on Matrix Analysis and Applications, 31 (2010) , pp. 1642-1662.

[0320]

[32] C.C. PAIGE AND M.A. SAUNDERS, Solution of sparse indefinite systems of linear equations, SIAM Journal on Numerical Analysis, 12 (1975) , pp. 617-629.

[0321]

[33] B.N. PARLETT, The Symmetric Eigenvalue Problem, Society for Industrial and Applied Math-ematics, 1998.

[0322]

[34] Y. SAAD, ILUM: a multi-elimination ILU preconditioner for general sparse matrices, SIAM Journal on Scientific Computing, 17 (1996) , pp. 830-847.

[0323]

[35] Y. SAAD, M. YEUNG, J. ERHEL, AND F. GUYOMARC'H, A deflated version of the conjugate gradient algorithm, SIAM Journal on Scientific Computing, 21 (2000) , pp. 1909-1926.

[0324]

[36] A. STATHOPOULOS AND J.R. MCCOMBS, PRIMME: Preconditioned Iterative Multimethod Eigensolver: Methods and software description, ACM Trans. Math. Softw., 37 (2010) , pp. 21: 1-21: 30.

[0325]

[37] A. STATHOPOULOS AND K. ORGINOS, Computing and deflating eigenvalues while solving multiple right-hand side linear systems with an application to quantum chromodynamics, SIAM Journal on Scientific Computing, 32 (2010) , pp. 439-462.

[0326]

[38] Y. SU, T. LU, AND Z. BAI, 2D eigenvalue problems I: Existence and number of solutions, arXiv preprint arXiv: 1911.08109, (2019).

[0327]

[39] U. VON LUXBURG, A tutorial on spectral clustering, Statistics and computing, 17 (2007) , pp. 395-416.

[0328]

[40] C. YANG, W. GAO, Z. BAI, X.S. LI, L. -Q. LEE, P. HUSBANDS, AND E. NG, An algebraic substructuring method for large-scale eigenvalue calculation, SIAM Journal on Scientific Computing, 27 (2005) , pp. 873-892.

[0329]

[41] H. ZHANG, B. SMITH, M. STERNBERG, AND P. ZAPOL, SIPs: Shift-and-invert parallel spec-tral transformations, ACM Transactions on Mathematical Software (TOMS) , 33 (2007) , pp. 111-127.

Claims

1. 1. A method for solving eigenvalues ​​of a matrix system utilizing a computing device, comprising using a processor of the computing device: loading the matrix system into a memory of the computing device; performing a partitioning into non-overlapping sub-regions and a rearrangement of the matrices of the matrix system; creating a boundary variable projection subspace associated with boundary variables of an adjacency graph of the matrix system, where the boundary variables are related to nodes of the adjacency graph that are connected to nodes located in neighboring partitions; creating an internal variable projection subspace associated with internal variables of the adjacency graph of the matrix system, where the internal variables are related to nodes of the adjacency graph that are connected only to nodes located in the same partition; constructing a projection matrix based on the boundary variable projected subspace and the internal variable projected subspace; determining eigenvalues ​​that solve a Rayleigh-Ritz eigenvalue problem using the projection matrix; and Utilizing said determined eigenvalues ​​in a technical application system. A method comprising:

2. The method of claim 1 , wherein the technological application system is selected from the group consisting of a clustering analysis system and a frequency response analysis system.

3. 3. The method of claim 2, wherein the clustering analysis system is an artificial intelligence (AI) classification system and the frequency response analysis system is selected from the group consisting of an audio analysis system and a building structural analysis system.

4. 4. The method of claim 1, wherein creating a boundary variable projected subspace associated with the boundary variables utilizes a zeroth order truncation of the spectral Schur complement matrix to avoid computing derivatives of the spectral Schur complement matrix.

5. Creating the boundary variable projected subspace associated with the boundary variables comprises: expanding the analytical eigenvectors of the spectral Schur complement matrix through their Taylor series; and injecting the leading eigenvector derivatives into the boundary variable projected subspace.

5. The method of claim 4, comprising:

6. 6. The method of claim 1, wherein the matrix system is a large sparse symmetric matrix pencil (A,M) that includes a matrix A and a matrix M, the large sparse symmetric matrix pencil constituting a domain of the matrix system.

7. The matrices A and M are symmetric, The matrix M is a symmetric positive definite matrix. The method according to claim 6.

8. The method of claim 6 , wherein the step of performing the domain partitioning utilizes an algebraic graph partitioner on an adjacency graph associated with the matrix system.

9. determining an eigenvector associated with each eigenvalue to form an eigenpair; and Utilizing said eigenpairs in said technical application system. The method of claim 1 , further comprising:

10. 16. An apparatus for an engineering analysis system comprising: memory; and 13. A processor comprising: loading a matrix system into said memory; performing a partition into non-overlapping sub-regions and a matrix reordering of the matrix system; creating a boundary variable projection subspace associated with boundary variables of an adjacency graph of the matrix system, where the boundary variables are related to nodes of the adjacency graph that are connected to nodes located in neighboring partitions; creating an internal variable projection subspace associated with internal variables of the adjacency graph of the matrix system, where the internal variables are related to nodes of the adjacency graph that are connected only to nodes located in the same partition; constructing a projection matrix based on the boundary variable projected subspace and the internal variable projected subspace; Utilizing the projection matrix to determine eigenvalues ​​that solve a Rayleigh-Ritz eigenvalue problem; and Utilizing said determined eigenvalues ​​in a technical application system. An apparatus comprising a processor configured to:

11. The apparatus of claim 10 , wherein the technological application system is selected from the group consisting of a clustering analysis system and a frequency response analysis system.

12. 12. The apparatus of claim 11, wherein the clustering analysis system is an artificial intelligence (AI) classification system and the frequency response analysis system is selected from the group consisting of an audio analysis system and a building structural analysis system.

13. 13. The apparatus of claim 10, wherein creating a boundary variable projected subspace associated with the boundary variables utilizes a zeroth order truncation of a spectral Schur complement matrix and avoids computing derivatives of the spectral Schur complement matrix.

14. Creating a boundary variable projected subspace associated with the boundary variables includes: expanding the analytical eigenvectors of the spectral Schur complement matrix through their Taylor series; and injecting the leading eigenvector derivatives into said boundary variable projected subspace; 14. The apparatus of claim 13, further comprising:

15. 15. The apparatus of claim 10, wherein the matrix system is a large sparse symmetric matrix pencil (A,M) that includes a matrix A and a matrix M, the large sparse symmetric matrix pencil constituting a domain of the matrix system.

16. The matrices A and M are symmetric, The matrix M is a symmetric positive definite matrix.

16. The apparatus of claim 15.

17. The apparatus of claim 15 , wherein performing the domain partitioning utilizes an algebraic graph partitioner on an adjacency graph associated with the matrix system.

18. The processor, determining an eigenvector associated with each eigenvalue to form an eigenpair; and Utilizing said eigenpairs in said technical application system.

18. The apparatus of claim 10, further configured to:

19. 1. A computer program for a technology analysis system, the computer program comprising: The procedure for loading a matrix system into memory; performing a partition into non-overlapping sub-regions and a rearrangement of the matrices of the matrix system; creating a boundary variable projection subspace associated with boundary variables of an adjacency graph of said matrix system, where said boundary variables are related to nodes of said adjacency graph that are connected to nodes located in neighboring partitions; creating an internal variable projection subspace associated with internal variables of the adjacency graph of the matrix system, where the internal variables are related to nodes of the adjacency graph that are connected only to nodes located in the same partition; constructing a projection matrix based on the boundary variable projected subspace and the internal variable projected subspace; determining eigenvalues ​​that solve a Rayleigh-Ritz eigenvalue problem using the projection matrix; and Utilizing said determined eigenvalues ​​in a technical application system. A computer program for executing the following:

20. the matrix system is a large sparse symmetric matrix pencil (A,M) that includes a matrix A and a matrix M, the large sparse symmetric matrix pencil constituting a domain of the matrix system; Matrices A and M are symmetric; The matrix M is a symmetric positive definite matrix; and performing the region partitioning step by utilizing an algebraic graph partitioner on an adjacency graph associated with the matrix system; 20. A computer program product as claimed in claim 19.

Images