Signal matrix solving pseudo-inverse function design method based on GPU calculation
By designing a GPU-based pseudo-inverse function for solving signal matrices, the problems of memory management and interface compatibility in large-scale matrix inversion are solved, achieving efficient parallel computing, which is applicable to fields such as signal processing and image reconstruction.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- EAST CHINA INST OF COMPUTING TECH
- Filing Date
- 2025-12-16
- Publication Date
- 2026-04-21
AI Technical Summary
Existing technologies suffer from fragmented memory management, high repetition, poor interface compatibility, and a single matrix storage method when performing large-scale matrix inversion, resulting in low computational efficiency and difficulty in meeting real-time requirements.
A pseudo-inverse function for solving signal matrices based on GPU computing is designed. By binding cusolverHandle and cublasHandle handles, it utilizes CUDA stream parallel computing, supports row/column master order matrix conversion, optimizes memory allocation and access, and follows the VSIPL interface specification to achieve parallel computing of pseudo-inverses of large-scale complex matrices.
It improves computational efficiency, reduces memory fragmentation, enhances interface compatibility, supports row/column master-order matrix input and output, and is suitable for scenarios with high real-time requirements such as signal processing and image reconstruction.
Smart Images
Figure CN121901552A_ABST
Abstract
Description
Technical Field
[0001] This invention belongs to the field of signal processing technology, specifically relating to a method for designing pseudo-inverse functions for solving signal matrices based on GPU computing. Background Technology
[0002] In the field of signal processing, finding the pseudo-inverse of complex matrices is widely used. VSIPL is a CPU-based computation function library containing many efficient operators for signal processing, and its interface design and usage specifications are standardized.
[0003] However, traditional CPU-based computation functions suffer from high computational complexity when dealing with large-scale matrices (pseudo-inverse calculations typically involve SVD decomposition, with a time complexity of O(n log n)). The insufficient computing power of CPU parallel computing leads to low computing efficiency, making it difficult to meet real-time requirements. Furthermore, the VSIPL standard library does not have a direct interface for finding the pseudo-inverse of complex matrices, requiring developers to repeatedly call matrix decomposition and matrix multiplication.
[0004] Furthermore, the existing GPU-based cuBLAS function library does not provide an interface for calculating pseudo-inverses. Developers need to use a combination of the GPU's cuBLAS library and the cusSOLVER library to implement matrix inversion, and all row-major matrices need to be manually transposed to column-major order, such as... Figure 1 As shown.
[0005] In summary, conventional methods have the following technical drawbacks: 1) Memory management fragmentation: Traditional methods use distributed memory allocation, which leads to low memory access efficiency, easily causes memory fragmentation problems, and affects computing speed.
[0006] 2) High repetition rate and poor interface compatibility: Existing GPU libraries do not provide a direct pseudo-inverse calculation interface, requiring repeated calls to the matrix decomposition and matrix multiplication interfaces for each calculation. Furthermore, incompatibility with the VSIPL specification makes it difficult to port projects implemented in VSIPL to GPUs.
[0007] 3) Limited matrix storage method: The existing cuBLAS and cuSOLVERDn libraries do not provide matrix output in row-major order, and the transpose of the output matrix needs to be implemented manually. Summary of the Invention
[0008] The objective of this invention is: 1. While adhering to the VSIPL interface specification, implement GPU parallel computation of large-scale complex matrix pseudo-inverses, facilitating the migration of VSIPL's CPU computation functions to GPU operations; 2. Optimize the memory allocation, access, and release efficiency of SVD decomposition and matrix multiplication; 3. Supports automatic conversion of row / column master order matrices to meet VSIPL specification requirements.
[0009] To address the aforementioned technical problems, the present invention provides a method for designing pseudo-inverse functions for solving signal matrices based on GPU computation, comprising the following steps: Obtain the complex matrix of the received signal from an N x M dimensional array containing the desired signal, interference signal, and noise, and having M snapshots and N array elements. The formula is as follows:
[0010]
[0011] in, It is the desired signal, the steering vector. Describes the phase / amplitude response of a signal on an array. It is the incident angle of the desired signal. It is an interference signal from a different direction. It's noise; Output in array With expected signal The optimal weighted coefficient vector is found under the premise of minimizing the mean square error. w The formula is as follows:
[0012] in, It is the pseudo-inverse of a complex matrix; The `void vsip_cmpinv_f` function, which takes `cusolverHandle` and `cublasHandle` handles as view pointers, binds different streams based on the `cusolverHandle` and `cublasHandle` handles. It also calls the `handle` parameter in the `cublas` library to bind the handle and bind independent CUDA streams, thus solving for the optimal weighted coefficient vector. w Computational tasks are migrated to run on GPUs; When solving, the complex matrix of the received signal from the N x M dimensional array is used. The mode type determines whether the matrix needs to be transposed. If so, the kernel function row2col_major_transpose_complex, which runs on the GPU, is used to convert it into a column-major ordered matrix. Based on the U, VT, and S matrices used in the SVD process, calculate the memory size occupied by lwork and the space required by the U, VT, and S matrices, and then calculate the total memory requirement of the column-major matrix. Implement single large-block memory allocation by designing the CudaResource class. For the S matrix, for each singular value ,satisfy If it is greater than the adaptive threshold ε, then The corresponding position is set to its reciprocal. Otherwise, it is 0, thus generating the S inverse matrix. The transpose of U yields the UT inverse matrix, and the transpose of VT yields the V inverse matrix. The pseudo-inverse of the complex matrix is obtained from the S inverse matrix, the UT inverse matrix, and the V inverse matrix. ; Output the pseudoinverse of the complex matrix in row-major or column-major order according to user requirements. Then, the optimal weighted coefficient vector w is obtained.
[0013] Preferably, based on the minimum mean square error criterion, the array output is... With expected signal The mean square error is the smallest.
[0014] Preferably, the array output With expected signal The formula for the mean square error is as follows: .
[0015] Preferably, the adaptive threshold ε is set to .
[0016] This invention provides a method for designing pseudo-inverse functions of signal matrices based on GPU computing. It follows the Vector Signal Processing Library (VSIPL) interface specification and supports both row-major and column-major matrix inputs. By optimizing singular value decomposition (SVD) and matrix multiplication operations, it achieves efficient parallel processing and computation of large-scale complex matrix pseudo-inverses of received signals on the GPU. This fills the gap in complex matrix inversion interfaces in the cuBLAS and VSIPL domains and is suitable for scenarios with high real-time requirements such as signal processing and image reconstruction. Attached Figure Description
[0017] Figure 1 A flowchart illustrating the steps for inverting an existing technology matrix; Figure 2 This invention provides a schematic flowchart of a method for designing pseudo-inverse functions for solving signal matrices based on GPU computing, as provided in an embodiment of the invention. Figure 3 This is a schematic diagram of multi-stream parallelism provided in an embodiment of the present invention; Figure 4This is a schematic diagram of unified memory management provided for an embodiment of the present invention. Detailed Implementation
[0018] The present invention will be further illustrated below with reference to specific embodiments. It should be understood that these embodiments are for illustrative purposes only and are not intended to limit the scope of the invention. Furthermore, it should be understood that after reading the teachings of this invention, those skilled in the art can make various alterations or modifications to the invention, and these equivalent forms also fall within the scope defined by the appended claims.
[0019] like Figure 2 As shown, this embodiment of the invention provides a method for designing a pseudo-inverse (Moore-Penrose inverse) function for solving a signal matrix based on a graphics processing unit (GPU), comprising the following steps: In the field of radar signal processing, the problem of inverting signal matrices is frequently encountered. For example, in adaptive beamforming, suppose there is a sensor array consisting of N elements that receives signals from M snapshots (i.e., sampled data at different times). Given the presence of the desired signal, interference signals, and noise, the mathematical model of the array's received signal can be expressed in the following matrix form: set up x Let (n) represent the signal vector (N x 1 dimensional complex vector) received by the N array elements at the nth snapshot. Then:
[0020] in, It is the desired signal, the steering vector. Describes the phase / amplitude response of a signal on an array. It is the incident angle of the desired signal; These are I interference signals from different directions; It's noise. For M snapshots, the received signals are combined into a matrix. X It is an N x M dimensional complex matrix, and its column vectors are the received signal vectors corresponding to each snapshot, i.e.:
[0021] Adaptive beamforming can be based on the minimum mean square error (MMSE) criterion, aiming to find a weighted coefficient vector. w (N x 1-dimensional complex vectors), such that the array output With expected signal The mean square error is minimized, that is:
[0022] By minimizing To find the optimalw We hope as close as possible From the perspective of the desired signal, for M snapshots, it can be written as:
[0023] in, Therefore, a solution is required. w We need to use the pseudo-inverse of a complex matrix, that is:
[0024] The function designed in this invention is named ( , , ,cublasHandle_t cublasHandle), where parameters A and A_pinv are the input and output matrices, and cusolverHandle and cublasHandle are the passed-in handles, which can be reused on the host code after creation. Figure 3 handle1, handle2, and handle3 are required parameters for calling the cublas library. Stream is a GPU-specific parallel computing stream module. By binding the stream and handle, the three matrices can be computed in parallel for SVD decomposition.
[0025] This function's interface design follows the VSIPL specification, and the VSIPL-CPU version code can be quickly migrated to GPUs for execution. The matrices to be computed are pre-allocated with device memory during host initialization to avoid impacting computational efficiency during memory allocation and matrix initialization. The matrix interface design is consistent with the VSIPL interface design specification, and the interface calls comply with VSIPL usage guidelines, facilitating the rapid migration of VSIPL-based computational tasks to GPUs. While standard interfaces in common VSIPL function libraries only require two view pointers for input and output matrices according to the VSIPL specification, this invention's pseudo-inverse computation function interface introduces new handle parameters of type `cusolverDnHandle_t` and `cublasHandle_t`. These handles are placed separately in the call parameters for convenient subsequent parallel computation and reuse. When performing inversions on a large number of matrices, creating handles internally within the interface has an unacceptable impact on computational efficiency. Therefore, this function design encourages the creation of multi-threaded matching CUDA streams and handles in the host code, binding an independent CUDA stream and handle (cublasHandle_t / cusolverDnHandle_t) to each matrix, supporting concurrent computation of multiple matrices. In the pseudo-inverse interface design, the added cublasHandle_t / cusolverDnHandle_t handles can be bound to different streams. When calling the pseudo-inverse interface, if multiple matrices need to be inverted, different streams can be set for each matrix. Streams can be reused, and similarly, cublasHandle_t / cusolverDnHandle_t handles can also be reused, so they can be used multiple times after being created in the main function. If cublasHandle_t / cusolverDnHandle_t is called within the interface, the handle needs to be recreated. Repeated creation of multiple handles leads to reduced function execution efficiency and consumes memory resources.
[0026] In the function implementation of void vsip_cmpinv_f, the kernel function row2col_major_transpose_complex is designed to automatically convert the row / column major order data of the input and output matrices in the parameters received by the function.
[0027] The vsip_cmpinv_f method of this invention can be used to solve problems quickly on a GPU. Thus, the solution is obtained. wIn the main program, the user creates stream1 and two handles: cublasHandle_t handle1 and cusolverDnHandle_t handle2. When vsip_cmpinv_f is called, the input matrix is X, and the output matrix is... Handle1 is named handle1, and handle2 is named handle2. The input is a complex matrix. X Following the VSIPL matrix standard specification, firstly, based on the matrix... X Initialize the specified mode type. The main order processing module determines whether the matrix needs to be transposed. When the mode of the input matrix is 0, the input matrix will be converted to column-major order for calculation. If the mode of the output matrix is 0, the matrix will be converted to row-major order for output.
[0028] The interfaces in the cublas library are all for column-major matrices and do not directly provide transpose functions. With this solution, users only need to provide the input and output matrix data in the interface parameters, without needing to worry about the internal implementation. Both the cusolver and cublas libraries output matrices in column-major order; therefore, a matrix order conversion module is added for cases requiring row-major output. Callers do not need to concern themselves with the matrix storage format; they only need to set the required storage format of the output matrix according to the VSIPL specification. This solution uses custom kernel functions "col2row_major_transpose" and "row2col_major_transpose" to implement matrix row and column transformations on the GPU, effectively utilizing the GPU's parallel computing capabilities.
[0029] The unified memory processing module calculates the required workspace and memory for the U, VT, and S variables generated by the economical SVD based on the input matrix, and then requests a unified memory pool total_memory.
[0030] The function implements single large-block memory allocation through the design of the CudaResource class. Using the SVD method to find the pseudo-inverse of a matrix involves the allocation and release of multiple small memory blocks. To avoid fragmentation caused by multiple memory allocations, this invention adopts a method of first calculating the total memory requirement and then uniformly allocating memory. Specifically, the required memory size is calculated based on the input matrix using the U, VT, and S matrices needed in the SVD process. The device memory size for parameters in the SVD calculation is also calculated, including the size of all device memory such as the diagonal matrix needed for subsequent matrix multiplication calculations. Memory alignment is performed on all calculated memory to optimize device memory access efficiency. The total memory requirement is allocated to a contiguous memory pool, and memory is used on demand. Furthermore, since GPU memory is used, resources are limited, and GPU data is reset if errors occur during GPU memory usage. Allocated memory must be initialized and released promptly. Therefore, this invention uses the RAII resource management class to uniformly initialize, allocate, and release memory, eliminating the need for manual management of each memory allocation and release, significantly reducing code complexity and improving code reusability. The execution steps of unified memory processing are as follows: Figure 3 As shown.
[0031] The CudaResource class supports Figure 4 All memory used during the computation is managed uniformly. `cusolverDnCgesvd_bufferSize` calculates the memory size of the lwork (temporary workspace) used during SVD computation, calculates the space required for matrices U, S, and VT, and sums them. A large memory pool is then allocated based on this total memory size. All subsequent memory is requested from this allocation pool, avoiding memory fragmentation. After all computations are complete, the allocated GPU memory does not need to be released by the user; the `CudaResource` class automatically releases the entire memory pool. Figure 4 The memory pool structure on the right shows the memory size occupied by each variable, as well as the starting address of the memory, base_ptr.
[0032] The adaptive SVD module employs an economical SVD approach to obtain the U, VT, and S matrices. Solving... The adaptive threshold is set to For each singular value ,satisfy If it is greater than ε, then The corresponding position is set to its reciprocal. Otherwise, it is 0.
[0033]
[0034]
[0035] The obtained V, After multiplying the UT matrices, we get .
[0036] =V× ×UT The main sequence conversion module outputs... The system will first check the data type and then output the data in either row-major or column-major order based on the user's requirements. Therefore, this invention can be used to migrate signal processing of adaptive beamforming from the CPU to the GPU, significantly improving computational efficiency.
[0037] This function uses economical SVD when calling the SVD interface, calculating only the singular vectors corresponding to non-zero singular values. This avoids calculating and storing useless null space vectors, significantly reducing memory usage, especially for large-scale non-square matrices. Singular value thresholding ensures numerical stability, allowing for a balance between accuracy and computational efficiency. The cusolverDnCgesvd interface in the cusolver library supports different types of singular vector calculation methods. This solution adopts the 'S' method provided by this function interface, called "economical SVD," which calculates only the left singular vectors of the first min (m,n) columns and the right singular vectors of the first min (m,n) rows. Economical SVD is mathematically equivalent to full SVD in pseudo-inverse calculations, but saves memory usage; for example, for a 512 x 400 matrix, it can reduce GPU memory resources by at least 20%. Furthermore, this invention uses the Thrust library to directly calculate the maximum singular value on the GPU based on a dynamic thresholding algorithm for determining the maximum singular value. This not only fully utilizes the parallelism of the GPU but also allows users to adjust the calculation precision, ensuring numerical stability.
[0038] The beneficial effects of the embodiments of the present invention are as follows: Excellent interface compatibility: The vsip_cmpinv_f interface is designed in accordance with the VSIPL specification, supporting row-major / column-major matrix input and output, and can be seamlessly integrated into existing VSIPL application frameworks. The interface parameter definitions are consistent with standard VSIPL functions. Although two handle parameters have been added, handle creation and initialization are simple, allowing for CPU-to-GPU computational migration without excessive modification to upper-layer application code, maintaining complete consistency of algorithm logic.
[0039] Memory utilization optimization: A unified memory allocation strategy is adopted, which reduces memory fragmentation by allocating large blocks of contiguous memory at once (such as merging and allocating data areas of matrices A, U, S, VH, etc.). Compared with distributed allocation, the unified allocation strategy reduces memory allocation operations by 9, thereby reducing the scheduling overhead of the CUDA driver, especially in high-density matrix computation scenarios. Furthermore, adaptive SVD reduces the use of video memory resources and improves computational efficiency by calculating only useful singular vectors.
[0040] Enhanced numerical stability: The adaptive singular value thresholding strategy (dynamically calculating tolerance = max_S × 1e-6 based on the maximum singular value) retains more effective singular values while ensuring numerical stability. This effectively avoids computational distortion caused by excessive truncation.
[0041] Significantly enhanced parallelism: The pseudo-inverse interface was designed with the simultaneous calculation of pseudo-inverses for multiple matrices in mind. By creating and reusing handles for both `cusolver` and `cublas` in the main function, combined with CUDA's multi-stream mechanism, the parallelism of pseudo-inverse operations is greatly improved. Furthermore, the reusability of handles reduces the overhead of handle creation and release.
Claims
1. A method for designing pseudo-inverse functions of signal matrices based on GPU computing, characterized in that, Includes the following steps: Obtain the complex matrix of the received signal from an N x M dimensional array containing the desired signal, interference signal, and noise, and having M snapshots and N array elements. The formula is as follows: in, It is the desired signal, the steering vector. Describes the phase / amplitude response of a signal on an array. It is the incident angle of the desired signal. It is an interference signal from a different direction. It's noise; Output in array With expected signal The optimal weighted coefficient vector is found under the premise of minimizing the mean square error. w The formula is as follows: in, It is the pseudo-inverse of a complex matrix; The void vsip_cmpinv_f function, which takes cusolverHandle and cublasHandle as pass-in handles as view pointers, binds different streams according to the cusolverHandle and cublasHandle pass-in handles, calls the handle parameter in the cublas library to bind handles, and binds independent CUDA streams, thus migrating the task of solving the optimal weighted coefficient vector w to the GPU for execution. When solving, the complex matrix of the received signal from the N x M dimensional array is used. The mode type determines whether the matrix needs to be transposed. If so, the kernel function row2col_major_transpose_complex, which runs on the GPU, is used to convert it into a column-major ordered matrix. Based on the U, VT, and S matrices used in the SVD process, calculate the memory size occupied by lwork and the space required by the U, VT, and S matrices, and then calculate the total memory requirement of the column-major matrix. Implement single large-block memory allocation by designing the CudaResource class. For the S matrix, for each singular value ,satisfy If it is greater than the adaptive threshold ε, then The corresponding position is set to its reciprocal. Otherwise, it is 0, thus generating the S inverse matrix. The transpose of U yields the UT inverse matrix, and the transpose of VT yields the V inverse matrix. The pseudo-inverse of the complex matrix is obtained from the S inverse matrix, the UT inverse matrix, and the V inverse matrix. ; Output the pseudoinverse of the complex matrix in row-major or column-major order according to user requirements. Then, the optimal weighted coefficient vector is solved. w .
2. The method for designing pseudo-inverse functions of signal matrices based on GPU computing as described in claim 1, characterized in that, Based on the minimum mean square error criterion, the array output is... With expected signal The mean square error is the smallest.
3. The method for designing pseudo-inverse functions of signal matrices based on GPU computing as described in claim 1, characterized in that, The array output With expected signal The formula for the mean square error is as follows: 。 4. The method for designing pseudo-inverse functions of signal matrices based on GPU computing as described in claim 1, characterized in that, The adaptive threshold ε is set to .