1
MANOJAVAM: A Scalable, Unified FPGA Accelerator for Matrix Multiplication and Singular Value Decomposition in Principal Component Analysis
arXiv:2605.01514v1 [cs.AR] 2 May 2026
Srivaths Ramasubramanian (Student Member, IEEE), Anjali Devarajan (Student Member, IEEE), Kousthub P Kaivar (Student Member, IEEE), Vibha Shrestta (Student Member, IEEE), Shashank D (Student Member, IEEE), Sowmyarani C.N (Member, IEEE), Govinda Raju M and K.S Geetha (Senior Member, IEEE)
Abstract—Principal Component Analysis (PCA) is widely used for dimensionality reduction in hyperspectral imaging, genomics, and neurosciences. However, it suffers from computational bottlenecks in matrix multiplication and singular value decomposition (SVD). Prior PCA hardware accelerators either target only one of these stages, rely on High Level Synthesis (HLS) that limits microarchitectural optimizations or use fixed point datapaths with limited dataset scalability. There is a need for a unified PCA accelerator that is suitable for datasets of any input dimension. Hence, the proposed work presents MANOJAVAM, a scalable PCA accelerator fabric, unifying matrix multiplication and SVD in a single architecture. MANOJAVAM(T,S) comprises an S number of TxT TPU-style systolic arrays employing block streaming for high-throughput matrix multiplication. It further integrates a highly parallel Jacobian unit implementing the Jacobi method for SVD with pipelined CORDIC based rotations. A two tier cache hierarchy and mode-aware memory policies adapts to the distinct memory access patterns of covariance matrix and rotation computation. For demonstration, MANOJAVAM(4,8) is realized on a Xilinx Artix-7 FPGA, achieving a frequency of 200 MHz at 1.271W. MANOJAVAM(16,32) is realized on Xilinx Virtex-Ultrascale+ FPGA, achieving a frequency of 434 MHz at 16.957W. Benchmarking on real-world datasets reveals that MANOJAVAM(16,32) achieves up to a 22.75x speedup in SVD latency and a 42.14x reduction in total energy consumption compared to a high-performance NVIDIA A6000 GPU. The architecture offers a unified, scalable, and energy-efficient platform for large-scale data analytics in both high-performance and edgecomputing environments. Index Terms—Hardware Acceleration, Domain Specific Architectures, Field Programmable Gate Arrays (FPGA), Principal Component Analysis, Singular Value Decomposition, Matrix Multiplication, Jacobi Algorithm, Energy Efficient Computing
I. I NTRODUCTION
T
HE end of Dennard scaling [1] and the slowing of Moore’s law [2] have made it difficult to meet the growing demands of modern day computational workloads. Modern day CPUs, though extremely fast in operation, are infeasible solutions to accelerate compute-intensive tasks due to This work was conducted at RV College of Engineering, Bangalore, India. Srivaths Ramasubramanian, Anjali Devarajan, Vibha Shrestta, Kousthub P Kaivar, Shashank D, Sowmyarani C. N, Govinda Raju M and K.S Geetha are with RV College of Engineering, Bangalore, India (e-mail: [email protected]; [email protected]; [email protected]; [email protected]; [email protected]; [email protected]; [email protected]; [email protected]).
their sequential nature. A major challenge in general purpose processors is that they are power hungry, as they offer support for multiple operations, have excessive precision support and a higher frequency of operations [3]. On the other hand, acceleration on Graphic Processing Units (GPUs) is largely constrained by very high power consumption and memory bottlenecks [4], thus limiting their full utilization of compute cores [5]. Thus, there is a need to accelerate such workloads with domain specific architectures [6]. As a result, hardware acceleration has emerged as a promising solution to speed-up massively parallel tasks and extract very high performance [5], [7]. Field Programmable Gate Arrays (FPGAs) have emerged as a primary platform for such tasks, offering gate-level reconfigurability and a high degree of parallelism. This enables rapid prototyping and the ability to implement highly optimized, domain-specific datapaths that are not constrained by the fixed instruction sets of general-purpose processors. Consequently, workloads accelerated on FPGAs show superior performance and energy efficiency compared to their general-purpose counterparts. Principal Component Analysis (PCA), is an important data mining technique used in several real world applications, such as image compression [8], image classification [9], facial detection [10], neurosciences [11] and genomics [12], [13]. PCA reduces the dimensionality of an input dataset by mapping it to a lower dimensional subspace. It involves the transformation of the dataset’s original features into a set of Principal Components, that maximally contain the variance of the input dataset while reducing the number of attributes [14]–[16]. The algorithm suffers from two bottlenecks: the computation of the Covariance matrix and Singular Value Decomposition (SVD). While software-based PCA implementations suffice for offline data analysis, they introduce non-deterministic latencies that are prohibitive for high-speed edge applications, such as autonomous drone navigation and real-time medical imaging. Existing general-purpose accelerators, such as mobile-GPUs, offer high peak throughput but often suffer from significant power overhead and memory bottlenecks during the iterative stages of eigendecomposition. Numerous attempts have been made to accelerate Principal Component Analysis [17]–[27], employing different methods to expedite the computation of covariance matrix and singular value decomposition. Prior work however, have suffered high resource consumption and operate on fixed input matrix di-
2
mensions. In this work, we present a parametric architecture where accelerator performance can be tuned by adjusting parallelism (S) and tile size (T ). While power dissipation naturally scales with these parameters, our baseline configuration of T = 4 and S = 8 achieves an low power profile of 1.271W on an Artix-7 FPGA, representing the lowest power consumption among the state of the art. The proposed work presents MANOJAVAM, a scalable and unified accelerator that integrates both matrix multiplication and SVD computation within a single architecture. The contributions of this work are as follows: • Unified Datapath: An accelerator that reuses its T × T systolic array cores for covariance matrix computation and Jacobi rotations, reducing resource utilization on FPGAs • Integrated Jacobi Unit with DLE: The Jacobi unit features a Data Lookup Engine (DLE) which directly interfaces with the outputs of the covariance matrix to obtain the pivot elements in a single scan once the accumulation is complete. This removes the need for round trips to BRAM to extract these elements. The DLE also implements tile-aware filtering to ensure it looks up at valid output tiles in the search of the pivot elements. • Numerical Robustness and Deterministic Convergence: A systematic Frobenius-norm-based evaluation across diverse datasets to establish a deterministic 50-sweep iteration schedule. This provides a significant safety margin for ill-conditioned datasets while eliminating the hardware overhead of complex on-chip convergence monitoring logic. • Mode-Aware Memory Subsystem: A specialized twotier cache hierarchy featuring two different cache writemiss policies to accommodate for distinct data patterns observed in covariance matrix computations and Jacobi rotations. • Multi-Platform FPGA Validation and Benchmarking: A comprehensive evaluation across the Artix-7 and the Virtex UltraScale+ FPGAs. Compared to a workstationgrade NVIDIA A6000 GPU, MANOJAVAM achieves up to a 28.2× speedup and over six orders of magnitude (1.3 × 106 ×) higher energy efficiency, while eliminating the non-deterministic software orchestration latencies native to GPU-based acceleration. II. L ITERATURE S URVEY Principal Component Analysis (PCA) remains a computational bottleneck in real-time and resource-constrained environments due to the high arithmetic and memory demands of covariance computation and singular value decomposition (SVD). A wide body of research has focused on accelerating portions of this pipeline using FPGAs, yet most designs treat the covariance and decomposition stages in isolation. This section traces the evolution of PCA hardware acceleration by progressively analyzing the architectural approaches to (A) covariance matrix computation, (B) SVD and eigendecomposition, and (C) unified designs that attempt to integrate both stages. We also survey (D) algorithmic strategies that
impact hardware parallelism and convergence, leading to (E) the motivation for MANOJAVAM. A. Covariance Matrix Acceleration The computation of the covariance matrix C = X T X, has drawn early attention in FPGA-based PCA accelerators. These designs often adopt block or vectorized matrix multiplication schemes to improve throughput and memory reuse. For instance, Korat et al. [17] presented an FPGA-based PCA pipeline that performs both learning and inference, using vector-matrix multiplication and QR decomposition for eigenspace extraction. While the design shows promise for small-scale matrices, its fixed-size datapaths and limited reuse strategy hinder scalability. Similar concerns are found in Mansoori et al. [19], who developed an HLS-driven architecture with dedicated Covariance and SVD units for hyperspectral imaging; while modular, its resource use grows steeply with dimensionality, and the high-level synthesis approach offers limited microarchitectural flexibility. A reconfigurable PCA accelerator was proposed by Shahrouzi et al. [21], which was capable of dynamically switching between matrix multiplication and QR stages. This dynamic reconfiguration was carried out by programming the FPGA with pre-compiled bitstream files for matrix multiplication and eigendecomposition. However, the need to reload bitstreams introduced significant reconfiguration overhead, resulting in performance degradation. Fernandez et al. [18] sidestepped this issue by offloading matrix multiplication to a host processor, integrating only the eigendecomposition in hardware. While this reduces hardware complexity, it prevents full integration of compute stages. Consequently, intermediate results must be transferred between configurations, leading to I/O overhead and reduced overall performance. Das et al. [20], on the other hand, applied PCA in FPGA-based network intrusion detection, where matrix multiplication was used to compute projections for dimensionality reduction. While their system demonstrates application-level effectiveness, it lacks on-chip decomposition logic and is tailored for fixed data characteristics, limiting its generality across domains. B. SVD and Eigendecomposition Acceleration Several works focus exclusively on accelerating the decomposition stage, often assuming a precomputed covariance matrix. These designs emphasize numerical stability, convergence speed, and hardware efficiency. Ma and Wang [28] leveraged the Hestenes-Jacobi method using fixed-point arithmetic and CORDIC engines, demonstrating fast convergence for matrices under 128×128. Yet the design’s reliance on external tools and point-wise rotation logic hampers scalability. Chen et al. [29] tackled the decomposition of dynamically varying matrices in wireless communication by designing a reconfigurable SVD ASIC for MIMOOFDM, integrating orthogonal reconstruction units. However, its serial rotation pipeline imposes latency constraints at scale. Wang and Zambreno [30] pursued a general-purpose FPGA SVD engine using IEEE-754 floating-point units and Hestenes-Jacobi iterations. While it supports large matrices
3
(up to 1024×1024), it lacks a matrix multiplication engine and relies on Xilinx IP blocks, making it less portable and energy efficient.
C. Unified Architectures for Covariance + SVD Despite the natural interdependence of the covariance and decomposition stages in PCA, few hardware accelerators unify both phases within a single architecture. The integration of matrix multiplication, decomposition, and projection within shared resources introduces complex synchronization, memory pressure, and latency management challenges. Among early efforts, Korat et al. [17] provided an end-toend FPGA design for small matrices, but with minimal parallelism or scalability. Mansoori et al. [19] connected Covariance and SVD units in a modular HLS design but did not optimize data movement or exploit systolic structures. Shahrouzi et al. [21] enabled time-sharing of logic between stages, but this led to pipeline stalls and suboptimal throughput. More structured integration appears in Zhang et al. [31], who used Brent-Luk systolic arrays and Jacobi rotations to construct a tightly coupled SVD engine on FPGA. While effective for 8×8 matrices, their design peaked at 98% resource utilization and could not scale to higher dimensions. Athi et al. [32] proposed a fast-converging parallel Jacobi method that accelerates convergence through strategic rotation ordering, yet the lack of a dedicated matrix multiplication core restricts its use in full PCA pipelines.
D. Algorithmic Strategies and Application-Level Contexts Beyond architectural integration, the choice of decomposition algorithm deeply influences convergence speed and hardware parallelism. Jacobi methods, while iterative, offer excellent orthogonality preservation and are naturally parallelizable — especially when applied via sweep-based scheduling. Torun et al. [33] systematically benchmarked Jacobi algorithms across CPU, GPU, and FPGA platforms. Their findings revealed significant memory irregularity on GPUs, resulting in poor cache utilization and warp divergence. FPGA-based designs, in contrast, performed well due to their predictable access patterns and deep pipelining capability. Athi et al. [32] further refined the Jacobi method by selecting maximal off-diagonal entries per sweep, achieving faster convergence and improved numerical stability — an approach that directly informs the rotation logic of our proposed design. Beyond architectural innovations, PCA has also been widely applied in various domain-specific systems as a pre-processing step. For example, PCA has been used for odor classification with neural networks [26], gas sensor signal processing [27], and data compression in wireless sensor networks [27]. Additionally, Ali et al. [27] proposed a hardware PCA implementation for gas identification using high-level synthesis on a Zynq SoC platform. While these works demonstrate PCA’s utility across domains, they are either application-driven or use high-level abstractions, and thus do not contribute to the architectural questions we address in this paper.
E. Summary and Motivation While a diverse set of designs target individual stages of PCA acceleration, there remains a lack of a cohesive architecture that combines high-throughput matrix multiplication, pipelined decomposition, and memory-aware data reuse in a fully scalable form. Existing solutions either operate on small matrices, omit key PCA stages, or depend on high-level synthesis and third-party IPs that limit control and scalability. To address these limitations, we propose MANOJAVAM — a unified, fully RTL-programmed architecture that tightly integrates systolic matrix multiplication with a pipelined Jacobi SVD engine. Our design introduces block streaming, operandmode-aware caching, and RTL-based sweep scheduling for parallel CORDIC-based rotations, all within a reconfigurable fabric. This combination of modularity, performance, and scalability sets MANOJAVAM apart from all prior works. III. P RINCIPAL C OMPONENT A NALYSIS Principal Component Analysis (PCA) is a method for dimensionality reduction that maps high-dimensional data onto a low-dimensional subspace, retaining maximum variance. It proceeds via the identification of orthogonal Principal Components (PCs), each formed as a linear combination of the original correlated features [14], [15]. PCA is described in Algorithm.1. Algorithm 1 Principal Component Analysis (PCA) Require: Dataset X ∈ RM ×N Ensure: Projected data O ∈ RM ×K 1: Standardize X 2: Compute covariance matrix C ← X ⊤ X 3: Compute eigenvalues and eigenvectors of C 4: Select top K components using EVCR or Scree plot 5: Project: O ← XVK 6: return O The initial step in the PCA pipeline involves dataset standardization, as defined in (1), to ensure that features having large numerical ranges do not disproportionately bias output principal components [14]. In the MANOJAVAM architecture, the input dataset X is assumed to be pre-standardized to zero mean and unit variance. This is done so that MANOJAVAM can entirely dedicate its datapath to accelerating covariance matrix computation and SVD, which are intensive operations and have time complexities of O(n · d2 ) and O(d3 ) respectively. yi =
xi − µ σj
(1)
Post standardization, the covariance matrix is constructed using equation (2). While C is symmetric, MANOJAVAM is designed to compute the full N × N covariance matrix to avoid complex control logic associated with computing only the upper or lower triangular matrix and reflecting the same to generate the full matrix. C = XT X
(2)
4
Following covariance matrix computation, the eigendecomposition is performed to extract orthogonal principal components. MANOJAVAM implements the Cyclic Jacobi Method for SVD acceleration due to its high degree of parallelism and suitability for fixed point precision [34], [35]. The hardware realization of this stage involves a specialized pipelined Jacobian Unit. First, a Data Lookup Engine (DLE) streams in the output covariance matrix to obtain the maximum off diagonal element cpq and its corresponding diagonal elements cpp and cqq . In order to achieve this through streaming in the covariance matrix at runtime, the DLE uses a combination of linear scan with tile aware filtering to obtain the maximum off diagonal element and corresponding diagonal elements on the fly. This approach helps avoid BRAM reads and passes to obtain the same elements, degrading performance. These elements are then streamed into the CORDIC Engine, which computes the rotation angle θ and subsequent trigonometric transformations. By decomposing the d × d eigendecomposition problem into a sequence of local 2 × 2 Givens rotations, the architecture maintains high numerical stability in fixed-point arithmetic. Upon the convergence of the Jacobi sweeps, the diagonal elements of the transformed covariance matrix represent the eigenvalues λ, which quantify the variance captured by each corresponding Eigenvector. To achieve dimensionality reduction, only a subset of these components are retained. The optimal number of retained components, denoted as k, is determined by analyzing the relative contribution of each Eigenvalue to the total variance. The choice of k can be determined easily using a combination of the Explained Variance Contribution Ratio (EVCR) and Cumulative Variance Contribution Ratio (CVCR). EVCR measures the individual contribution of the ith Eigenvector to the pool of Eigenvectors, and CVCR quantifies the total variance retained when mapping the data into the space encompassed by the first k principal components. The equation for EVCR and CVCR is given in (3) and (4). EV CR = PN
λi
m=1 λm
Pk λi CV CR = PNi=1 m=1 λm
(3)
(4)
Where λi denotes the ith Eigenvalue obtained. Another method frequently employed is the Scree plot, in which the Eigenvalue index is plotted along the X-axis, and the Eigenvalue along the Y-axis. The plot is always a decreasing plot. The rule of thumb is to usually consider one less of all the Eigenvalues to the left of the inflection point in the curve, although this must be combined with EVCR or CVCR for robust analysis. Finally, the original dataset is then projected onto the new subspace. This is achieved after the eigenvectors are arranged in columns, and the original dataset matrix is multiplied by the Eigenvector matrix. (5) depicts the projection. Om×k = Xcentered Vk ∈ RM ×K
(5)
IV. C OMPUTATIONAL B OTTLENECKS A SSOCIATED WITH PCA The primary challenge in scaling Principal Component Analysis (PCA) for high-dimensional, real-time datasets lies in the asymmetric growth of its core computational stages. While the covariance matrix computation scales at O(n · d2 ), the subsequent eigendecomposition (via the Cyclic Jacobi method) requires O(d3 ) operations per sweep. Since the Jacobi method typically converges within a constant number of sweeps, the O(d3 ) complexity dominates only when the feature dimension d is large relative to the sample size n. To establish a high-performance baseline, the computational bottlenecks of the PCA pipeline were profiled on an NVIDIA RTX A6000 GPU (Ampere architecture, 48GB GDDR6 VRAM). As illustrated in Fig. 1a and Fig. 1b, the execution time was analyzed across two distinct scaling regimes. In the constant row regime (Fig. 1a), it is observed that as the number of features in the dataset d approaches 1000, the O(d3 ) complexity of the Jacobi rotations begin to dominate the total latency. Conversely, in the constant feature regime (Fig. 1b), increasing the datapoints n to 105 shifts the bottleneck towards covariance matrix computation. These observations necessitate the need for an accelerator to accelerate covariance matrix computation and SVD for PCA. Unlike GPUs that suffer from non-linear scaling and software-stack overheads, our unified FPGA design provides a deterministic, pipelined dataflow that accelerates both matrix multiplication and iterative rotations in a single, highthroughput fabric. V. T HE JACOBIAN A LGORITHM The Jacobi algorithm is widely employed to compute the eigendecomposition of the covariance matrix. When compared to the Golub-Kahan algorithm [36], [37], the Jacobi algorithm presents the following advantages: It is highly parallelizable, results in symmetric memory access patterns, and supports fixed point implementations [34]. Furthermore, it avoids the hardware-intensive divisions and square roots required by the QR methods, replacing them with area-efficient CORDIC micro operations [38], [39]. The Jacobi algorithm is presented in Algorithm.2. The algorithm proceeds by obtaining the largest off-diagonal element cpq of the covariance matrix C and the corresponding indices p and q. The rotation angle, θ, is then calculated employing (6). 2cpq 1 −1 (6) θ = tan 2 cpp − cqq The rotation angle θ reorients the coordinate system such that element cpq becomes zero, resulting in the covariance matrix moving closer to its diagonal form. The Data Lookup Engine (DLE) executes a pivoting strategy by identifying the maximum off-diagonal element (ai,j ) for each rotation. While the CORDIC units wait for these indices, this approach ensures that each iteration achieves the maximum reduction in offdiagonal energy.
5
(a)
(b)
Fig. 1. Breakdown of PCA execution time into matrix multiplication and SVD components under different dataset dimensions. (a) SVD dominates when the number of features increases. (b) Matrix multiplication dominates as the number of rows increases.
Algorithm 2 Jacobi Eigenvalue Algorithm Require: Symmetric matrix C ∈ RN ×N Ensure: Diagonalized matrix C, eigenvectors in V 1: Initialize V ← I 2: while off-diagonal norm > ε do 3: Find indices (p, q) suchthat |cpq| is maximum 2cpq 4: Compute θ ← 12 tan−1 cpp −c qq 5: Set cos θ and sin θ accordingly 6: Construct rotation matrix R as identity with: 7: Rpp ← cos θ, Rqq ← cos θ 8: Rpq ← sin θ, Rqp ← − sin θ 9: Update C ← R⊤ CR 10: Update V ← V R 11: end while 12: return C, V
Since the covariance matrix C is real and symmetric, the spectral theorem applies directly, enabling its diagonalization via orthogonal similarity transformations, depicted in (10). ′
C = RT CR
The Jacobi algorithm applies a sequence of such rotations—referred to as Jacobi sweeps—to iteratively zero out off-diagonal elements. The method exhibits quadratic convergence, whereby each sweep significantly reduces the magnitude of the largest off-diagonal entry. The convergence of the Cyclic Jacobi method is typically monitored using the off-diagonal Frobenius norm, defined as: v uX n u n X |aij |2 (11) Eof f (A) = t i=1 j̸=i
The rotation matrix R is arrived at by initializing an identity matrix, and populating Rpp = cosθ, Rpq = sinθ, Rqp = −sinθ and Rqq = cosθ as shown in (7). 1 0 ··· ··· 0 0 1 · · · · · · 0 .. .. . . .. . . . . cos θ sin θ R= (7) − sin θ cos θ . . .. .. .. .. . . 0 0 ··· ··· 1 The Jacobi algorithm is grounded in a similarity transformation framework, which expresses the diagonalization of a matrix A through a sequence of orthogonal transformations, as shown in (8). ′ A = P −1 AP (8) When the transformation matrix P is orthogonal (ie., P −1 = P ), this reduces to the Spectral theorem, as described in (9). T
′
A = P T AP
(10)
(9)
where aij represents the off-diagonal elements of the evolving covariance matrix. In many software-based SVD solvers, an ’early-exit’ strategy is employed where the algorithm terminates once Eof f falls below a predefined threshold ϵ. However, implementing a real-time Frobenius norm unit on-chip presents significant hardware challenges. The calculation requires a high-precision Square-Root-of-Sum-of-Squares (SRSS) pipeline across the entire N × N matrix, which introduces substantial logic overhead, routing congestion, and a potential reduction in the maximum clock frequency (Fmax ). To maintain the high-throughput, streaming nature of MANOJAVAM, we move the convergence analysis offline. By performing a comprehensive Frobenius norm study across varied data modalities, we determine a fixed iteration count (i) that guarantees convergence. The fixed iteration count is also made large enough to accommodate ill-conditioned data at the input of the accelerator. This approach eliminates the need for complex monitoring logic while ensuring a deterministic execution latency. VI. A RCHITECTURE OF M ANOJAVAM This section introduces MANOJAVAM, a hardware accelerator, to speed up matrix multiplication and SVD for PCA.
6
The core of MANOJAVAM’s architecture is the Matrix Multiplication Engine (MM-Engine), which integrates an array of T xT systolic arrays alongside dedicated matrix accumulators. This unit performs both covariance matrix computation and rotations by employing block streaming to scale to large input matrices. Eigendecomposition is enabled using a lightweight Jacobian Unit, which identifies the maximum off-diagonal entry of the covariance matrix, and employs CORDIC-based micro-operations to calculate the parameters of the Given’s rotation matrix. Rotations are carried out in the MM-Engine, eliminating the need for a separate rotation datapath, and reducing hardware redundancy. A two-level cache hierarchy ensures optimized operand delivery, and a coordinated hierarchy of RTL controllers ensures synchronization across all datapath units. The high-level architecture of MANOJAVAM is shown in Fig.2. Central to its efficiency is a unified data-path that treats the covariance computation and eigendecomposition not as separate tasks, but as two operational modes of the same high-throughput fabric. A. Matrix Multiplication Accelerator Architecture MANOJAVAM’s matrix multiplication accelerator adopts a block streaming methodology. Partial product submatrices are computed in parallel and accumulated across corresponding row and column blocks to form the complete matrix product [19].This streaming strategy guides the architectural design of the matrix multiplication core, comprising an array of compact T × T systolic arrays. To mitigate the memory wall encountered by general-purpose GPUs (as demonstrated in Section IV), MANOJAVAM utilizes a deterministic blockstreaming methodology. By partitioning large N × D matrices into T × T tiles, the architecture maintains a constant memory footprint regardless of the sample count, ensuring that the system avoids the latency penalties associated with cache thrashing. Each systolic array in the MM-Engine is responsible for computing a submatrix of the final product matrix. Compared to a monolithic Tensor Processing Unit (TPU) core, this modular approach provides several key advantages: support for arbitrary matrix dimensions, efficient handling of boundary cases, improved resource utilization, and reduced power consumption via selective unit-level gating. In contrast, a large TPU becomes suboptimal when input matrix dimensions are misaligned with its array configuration, leading to under-utilization with small matrices and complex batching strategies for larger inputs [40]. The modular approach allows for efficient processing of rectangular matrices, commonly seen in PCA workloads. Block streaming is highly effective when dealing with very large matrix operands that cannot be streamed in all at once into the core computing units. A ”block” of tiles is streamed sequentially to compute the partial product tiles. These tiles are then accumulated to create the output tile [41]. An illustration of block streaming is shown in Fig.3. Systolic arrays are widely adopted for implementing matrix multiplication in hardware due to their highly regular, pipelined structure, enabling efficient spatial dataflow and minimal control overhead. Their inherent regularity allows for
compact layouts and predictable timing, making them wellsuited for FPGA and ASIC implementations. Systolic architectures have been employed in a variety of high-performance systems, most notably in Google’s Tensor Processing Unit (TPU) [42] and several other prior work on Machine Learning accelerators, DNN Accelerators and Digital Signal Processing [43]–[53]. Compared to traditional row-column multipliers, systolic arrays offer high parallelism, better resource utilization, and improved scalability for large matrix operations. As depicted in Fig.4, input operands are streamed into the systolic array in a systolic fashion, as suggested by the parallelogram input data profile to the systolic array. Each element of the systolic array is a Multiply and Accumulate (MAC) unit. Data is fed into these Processing Elements (PEs) in a skewed systolic fashion, allowing for 100% functional occupancy of the MAC units once the pipeline is primed. This spatial dataflow is architecturally significant as it minimizes global wire congestion and reduces power consumption by localizing data movement strictly between adjacent PEs, avoiding the high-energy cost of global bus transactions. Matrix multiplication is performed over multiple passes, with the number of passes determined by the dimensions of the input operand matrices. The MM-Engine in MANOJAVAM consists of S T xT systolic arrays, enabling up to S matrix multiplication operations to be parallelized concurrently. Each systolic array is responsible for computing a specific submatrix of the overall matrix product, corresponding to a pair of row and column blocks. These row and column blocks are further partitioned into T xT tiles that align with the dimensions of the systolic arrays. This tile-based partitioning supports highthroughput block streaming and ensures that data is processed in a spatially and temporally efficient manner. For instance, in computing the covariance matrix C, the submatrix C00 is calculated by the systolic array SA0 , which accumulates the products of row block-0 of X T and column block-0 of X. Each systolic array is paired with a dedicated accumulator to maintain correct accumulation of partial products across passes. This preserves numerical consistency and prevents accumulation of the incorrect matrix product from the MMEngine. First, the tiles are queried from the shared LHS and private RHS caches associated with the systolic array-matrix accumulator pair. The tiles first pass through Matrix Padding Units (MPU) at the interface between the MM-Engine and the caches. The MPU aids in inputting the tiles in a systolic fashion. There are two separate MPUs developed for the two input operands, as they have different systolic profiles. After the formation of matrix products, the output from the systolic array is fed into the matrix accumulator. This accumulates the T xT product tiles associated with the row and column block. Upon completion of all row and column tile iterations, the toplevel controller signals the end of partial product submatrix computation. The output from the matrix accumulator is then forwarded both to the Jacobi Unit for eigendecomposition and to memory for covariance matrix storage. The MM-Engine is reused for performing matrix multiplication during rotations. Crucially, the MM-Engine is not idle during the eigendecomposition phase; it is repurposed to apply the
7
Fig. 2. Manojavam : High Level Architecture
Fig. 3. Block Streaming Illustration
calculated Givens rotations to the entire covariance matrix. By acting as a parallel transformation engine that updates multiple rows and columns simultaneously, MANOJAVAM eliminates hardware redundancy and maximizes the area-efficiency of
Fig. 4. Systolic Array
8
the silicon footprint. To manage this distinction, the top-level controller issues a one-bit mode signal to the datapath, indicating whether the system is performing covariance computation or rotation at a given time. This ensures that the rest of the datapath executes the correct operation while appropriately scheduling operands for matrix multiplication during both covariance computation and rotation stages. This allows the reuse of computational resources of the accelerator across both stages. 1) Illustration: To illustrate the working of the matrix multiplication accelerator, assume an accelerator configuration of T = 4 and S = 8. Consider the input dataset to be of size 1000x1024, ie-1000 rows and 1024 features. If this is represented by matrix X, then X is of dimension 1000x1024, and X T is of dimension 1024x1000. Thus, the number of row blocks in X T is 1024/4 = 256, and the number of column blocks in X is again 1024/4 = 256. Each row and column block is divided into 1000/4 = 250 tiles, each of dimension 4x4. In pass 1, the 8 systolic arrays are scheduled as follows - SA0 is associated with the computation of partial product matrix R0 C0 , SA1 is associated with the computation of partial product matrix R0 C1 , SA2 is associated with the computation of partial product matrix R0 C2 and similarly, SA7 is associated with the computation of partial product matrix R0 C7 . This is held on until all 256 tiles associated with R0 , and each of the 256 tiles associated with C0 to C7 , are scheduled and their output product tiles are computed and forwarded to the matrix accumulators indexed with the systolic array. After all tiles in the row block and all of the column blocks have been scheduled, the top-level controller asserts that the 8 partial product submatrices have been successfully computed, and it arranges for the subsequent row blocks and column blocks to be fed in. As there are 256 column blocks to be computed, the top level controller schedules C8 to C15 for computation, while retaining R0 . The process continues until all of the 256 column blocks have been passed over for R0 . This ensures that one row of the output covariance matrix is completed, and continues until all the row blocks are iterated. B. Cache Subsystem MANOJAVAM employs a two-tier cache hierarchy system tailored for efficient matrix operand delivery. Operand A (LHS) is served by a single shared cache while operand B (RHS) is distributed across S private caches. Each of the private caches is tightly coupled to its corresponding systolic array and matrix accumulator. In the architecture, the shared and private caches are directly mapped. This cache system organization is architecturally driven - operand A is broadcasted and reused across multiple passes across systolic arrays, while operand B is unique to each array instance. These data access patterns have led to the shared and localized caching to ensure efficient operand delivery. This asymmetric cache architecture is specifically engineered to facilitate high efficiency block streaming. By utilizing a shared LHS cache, the architecture performs a single ’broadcast’ read, serving all processing elements simultaneously and reducing global memory bandwidth requirements by a factor of S. Conversely,
each systolic array is equipped with a Private RHS Cache to store the distinct vertical tiles required for its specific submatrix computation. Unlike conventional cache designs that store matrix rows directly, MANOJAVAM cache store complete matrix tiles, each encompassing a T × T submatrix in a single row entry. This layout avoids the inefficiency of issuing T separate cache reads to reconstruct one tile. Instead, a tile is fetched entirely in a single burst read to cache memory. An offline script flattens the raw input dataset into the proposed layout and then loaded into memory at system reset. The cache rows are populated upon cache misses. Each cache block is managed by a cache controller, totaling the number of controllers across the architecture to (S + 1). These cache controllers independently issue address instructions for operand data from memory. These controllers dynamically adjust their cache write-miss policies based on the phase of PCA computation. When computing the covariance matrix, C = X T X, the access pattern resembles the Livermore loop benchmark [54], given in (12). Hence, during the computation of covariance matrix, the caches must operate with the writearound policy for cache misses to ensure no cache pollution and high performance [55]. A[j] = B[j] + C[j]
(12) ′
T
In contrast, during the rotation phase, ie- C = R CR and V = V R, the access pattern aligns with the SAXPY-like behavior in the Linpack benchmark [56], given in (13), where write-allocate no-fetch-on-write policy is more appropriate [55]. A[j] = A[j] + B[j] (13) A mode signal propagates through the system to reconfigure cache controller behavior at runtime, enabling seamless adaptation between these two execution modes and ensuring policy coherence with workload characteristics. C. Jacobian Unit The Jacobian Unit in the accelerator serves as the bridge between the covariance matrix and the computation of the Givens rotation matrix, R. The Givens matrix is employed in rotations to converge at the eigenvalues and eigenvectors of the input covariance matrix. The high level design of the Jacobian unit is given in Fig. 5. The Jacobian unit comprises of the Data Lookup Engine (DLE), which interfaces directly with the output ports of the systolic arrays. The DLE streams in the covariance matrix directly from the accumulation stage and uses a linear scan to identify the maximum off-diagonal element cpq , corresponding indices (p, q) and elements cpp and cqq for the CORDIC engine. By resolving the pivot coordinates at the hardware interface, the DLE removes the need for redundant memory round-trips to the BRAM, providing a low-latency transition to the rotation phase. To maintain mathematical validity during the search, the Jacobian Controller implements tile-aware filtering. By monitoring the current row-block index, the DLE selectively masks
9
Fig. 5. Jacobian Unit Architecture
the main diagonal elements (cpp ) during the comparison phase, ensuring the engine only locks onto valid off-diagonal candidates. On reset, a global register is initialized to track the maximum off-diagonal element and its associated diagonal terms cpp and cqq . As each matrix multiplication pass produces S partial results from the accumulators, the Jacobian Controller inspects the current row block index and selectively filters out diagonal elements accordingly. For instance, during the processing of row block R0 , the diagonal elements from Acc0 (which correspond to the main diagonal of C) are discarded, while all elements from accumulators Acc1 –AccS are considered. In subsequent passes within R0 , all elements are valid for checking. A similar rule applies for row block R1 , where only Acc1 diagonals are excluded. This localized, tileaware filtering ensures that only valid off-diagonal candidates are passed to the DLE. Once all checks are completed for a given row block, the controller outputs the values of cpq , cpp and cqq , along with their corresponding indices p and q, for use in the Givens rotation step. The outputs from the DLE—specifically the values cpq , cpp and cqq , are forwarded to the CORDIC kernel [38], [39], which computes the sine and cosine of the rotation angle θ required for the Jacobi transformation. The rotation angle is computed using the relation in (6). This is implemented via a pipelined CORDIC arctangent unit, followed by a 1-bit right shifter. The resulting angle is then fed to additional CORDIC units that compute sin(θ) and cos(θ) in parallel. These trigonometric values are then passed to the Givens Controller, which initializes the corresponding entries of the Givens rotation matrix in memory. This matrix is pre-initialized to the identity matrix upon reset, and subsequently updated at runtime. Since the Givens matrix is reused in multiple operations, a write-back to memory is necessary and architecturally justified. The Givens
Controller uses the indices p and q to identify the appropriate memory locations, and updates them in a single event using a dual-port memory system. Once written, the top-level controller signals the beginning of the rotation phase by asserting the mode signal. Although initial accesses to the rotation matrix result in cache misses, the tiles are populated in the LHS and RHS operand caches. VII. E XPERIMENTAL S ETUP AND R ESULTS A. FPGA Implementation MANOJAVAM is developed entirely in synthesizable Verilog (IEEE Std 1364-2005) using the AMD Vivado 2024.1 design suite. The design is implemented on two platforms - a Xilinx Artix-7 xc35atcpg236-1 FPGA and Xilinx Virtex Ultrascale+ xcu250-figd2104-2L-e FPGA. The platforms are chosen to showcase design’s deployability across both resource-constrained and resource expansive targets. The ample resources in the Virtex Ultrascale+ FPGA help support the conduct of ablation studies for different tile sizes (T ) and parallelism indices (S). Functional correctness was validated using custom simulation testbenches, while Xilinx Design Constraints (XDC) were applied to guide placement and routing through floorplanning. Upon implementation, the following utilization report estimates were obtained, and are depicted in Table.I and Table.II. TABLE I M ANOJAVAM (4,8) FPGA R ESOURCE U TILIZATION - X ILINX A RTIX -7 XC 35 ATCPG 236-1 Resource LUT FF BRAM DSP
Usage 9796 23077 30.5 64
Total Available 20800 41600 50 90
Utilization 47.10 55.47 61 71.11
10
TABLE II M ANOJAVAM (16,32) FPGA R ESOURCE U TILIZATION - X ILINX V IRTEX U LTRASCALE + XCU 250- FIGD 2104-2L- E Resource LUT FF BRAM DSP
Usage 195814 143777 940.5 4096
Total Available 1728000 3456000 2688 12288
Utilization 11.33 4.16 34.99 33.33
MANOJAVAM operates at 200 MHz with 1.271 W power dissipation on the Artix-7 FPGA, and at 434 MHz consuming 16.957 W on the Virtex Ultrascale+ platform. To provide a conservative performance bound, we developed a cycle-approximate analytical simulator that models a worstcase sequential dataflow. This model accounts for effective access times (EAT) by incorporating a cache hit rate of p = 0.9 and a 10× penalty for off-chip DRAM access. For a tile size of T = 16, the simulator calculates the total execution time as the aggregate of data-loading overhead and systolic computation cycles. This ensures that the reported performance metrics represent a strictly attainable lower bound, even under significant memory contention. The Table.III outlines the comparison between prior accelerators and MANOJAVAM. B. Execution Time Analysis on Real World PCA Benchmarks This section evaluates the execution time performance of MANOJAVAM on a suite of real-world benchmark datasets commonly employed in Principal Component Analysis. The selected datasets summarized in Table.IV span diverse data modalities and dimension characteristics. These datasets cover a range of application domains including computer vision, hyperspectral imaging, biomedical analysis and text mining. For comparisons, MANOJAVAM is benchmarked against an NVIDIA A6000 GPU [65]. The total execution times for various datasets across all platforms is presented in Fig.6. The results demonstrate that MANOJAVAM outperforms the GPU on all datasets. Notably, on the CIFAR-10 dataset which is the most compute intensive benchmark in our dataset suite, MANOJAVAM achieves a 3.87× reduction in total execution time compared to the NVIDIA A6000. The benchmarking results from the MNIST-8x8 and Breast Cancer datasets show that the GPU performs sub-optimally compared to MANOJAVAM. This is attributed to the small sizes of the covariance matrix generated from these datasets which fail to efficiently utilize the massive parallel fabric of the NVIDIA A6000 [33]. The performance bottleneck on the GPU is driven by kernel launch latencies and SIMT (Single Instruction Multiple Thread) branch divergence during iterative Jacobi sweeps. In contrast, MANOJAVAM avoids these overheads through dedicated, hard-wired control logic and a parametric tiled dataflow, ensuring high utilization regardless of matrix scale.
consumption (E) is defined as the product of peak measured power (Ppeak ) and the total end to end execution latency (Ttotal ). The results show that MANOJAVAM showcases a significant reduction in energy consumption across all the datasets compared to the NVIDIA A6000. Notably, for smaller datasets such as the MNIST-8x8, the energy gain exceeds five orders of magnitude (> 105 ×). This is attributed to a higher power floor in the GPU and CUDA overhead times resulting in higher energy consumption, while the proposed accelerator has a much lower power floor and is purely datapath, resulting in low power operation. In high-dimensional scenarios such as CIFAR-10 (N = 3072) and 20 Newsgroups (N = 1024), where computational kernels typically saturate GPU resources, MANOJAVAM continues to exhibit superior efficiency. The MANOJAVAM(16,32) configuration, deployed on the Virtex UltraScale+, achieves a 42.14× energy reduction for CIFAR-10 compared to the A6000. This confirms that the proposed systolic Jacobi array and pipelined matrix multiplication units provide higher computational density per Watt than the SIMT-based kernels of modern workstation GPUs. The results further highlight the advantage of cycledeterministic dataflow. Unlike the A6000, which exhibits power fluctuations due to dynamic frequency scaling and nondeterministic OS-level interrupts, MANOJAVAM provides a stable, predictable power profile. This makes the proposed accelerator particularly suitable for deployment in missioncritical edge environments where energy budgets are stringent and thermal throttling must be avoided. D. Frobenius Norm Convergence Analysis To validate the choice of the number of Jacobi sweeps, we analyzed the relative Frobenius norm convergence across diverse datasets: MNIST, Olivetti Faces, and Breast Cancer. As illustrated in Fig.8, the relative off-diagonal energy for most standard datasets saturates at the numerical noise floor within 10 to 15 iterations.
C. Energy Efficiency and Computational Density Fig.7 presents a comparison of the energy profiles for the various datasets across the hardware platforms. The energy
Fig. 8. Relative error (Eof f ) vs. number of sweeps for multiple datasets
11
TABLE III C OMPARISON WITH P RIOR PCA ACCELERATORS Architecture LUT FF BRAM DSP Fmax (MHz) Power (W) Maximum Dimension (AxB) Korat.et.al [17] 272026 196 2464 183 16x30 Mansoori.et.al [19] 47880 34048 38 117 90 2.37 640x480 Shahrouzi.et.al [21] 24692 49384 11 100 1.58 64x3823 Wang.et.al [30] 13372 12723 129.2 8x8 Das.et.al [20] 17971 8932 74 139.5 Torun.et.al [33] 178 0.26 8x8 Fernandez.et.al [18] 16045 14417 28 124.6 302500x224 Ma.et.al [28] 182605 3044 740 183 128x128 Shiuping.et.al [31] 63656 63656 130 8x8 Athi.et.al [32] 32368 18760 48 192 236 8x8 Bravo.et.al [57] 40 43 112.4 256x256 Kasap.et.al [58] 23567 21268 53 64 103.88 4x4 MANOJAVAM (4,8) 9796 23077 30.5 64 200 1.271 Scale-Invariant MANOJAVAM (16,32) 195814 143777 940.5 4096 434 16.957 Scale-Invariant Note: “–” indicates that the corresponding figure was not reported in the referenced work. Note: Unlike prior works constrained by fixed on-chip buffer sizes, MANOJAVAM’s block-streaming allows processing of arbitrarily large matrices limited only by external storage capacity. TABLE IV S UMMARY OF B ENCHMARK PCA DATASETS Dataset MNIST-8x8 [59]
Number of Records 1797
Number of Features 64
MNIST-28x28 [60]
70000
784
CIFAR-10 [61]
60000
3072
Olivetti Faces [62]
400
4096
Breast Cancer [63]
45312
7
20-Newsgroups [64]
18846
1024
Despite this rapid convergence in typical cases, MANOJAVAM maintains a fixed limit of 50 iterations as a safety factor. This high upper ceiling is chosen to handle illconditioned datasets—where eigenvalues are closely clustered and require more rotations to achieve separation. This design choice ensures that the accelerator provides a universal ”Factor of Safety,” guaranteeing variance accuracy even for illconditioned input data distributions. Consequently, the architecture trades a negligible amount of redundant computation for a significant gain in hardware simplicity and cross-dataset reliability. VIII. D ESIGN S PACE E XPLORATION MANOJAVAM’s two tunable parameters include the Tile Size (T ) and the Parallelism Index (S). The selection of the two parameters has a drastic impact on the performance of the accelerator. While it is favorable to choose a very large value of S and T to extract the highest possible execution speed, it has negative implications on resource consumption, power dissipation, and size. Thus, it is crucial to choose S and T based on the application and the size of the input dataset. MANOJAVAM, being scalable, can be tailored to meet any PCA workload given its tunable tile size and cores deployed, making it extremely versatile.
Brief Description Grayscale images of handwritten digits (8×8 resolution), widely used for digit classification tasks. Standard MNIST dataset with 28×28 grayscale images of digits; high-dimensional image input. RGB images (32×32×3) across 10 object categories; common for low-resolution object recognition. Grayscale face images (64×64 pixels) of 40 individuals with varying facial expressions and lighting. Biomedical features extracted from mammographic scans for breast cancer diagnosis. Text documents represented as TF-IDF feature vectors across 20 categories; high-dimensional sparse input.
The following subsections describe the variation in performance of the proposed accelerator on execution time, power consumption and resource utilization across different tile sizes and cores deployed. These ablation studies are performed on the Xilinx Virtex Ultrascale+ xcu250-figd2104-2L-e FPGA. A. Execution Time Analysis for Varying Tile Sizes The execution time of the accelerator varies inversely with the square of the tile size (T 2 ). This is due to the fact that the input dataset of dimensions (M, N ) is split into tiles, resulting in MT N 2 tiles of computation. This is depicted in Fig.8(a), where the execution time is profiled for varying tile sizes across benchmark datasets for a fixed value of S = 32. B. Execution Time Analysis for Varying Parallelism Index The execution time of the accelerator varies inversely with the parallelism index (S). This is due to the fact that the input dataset of dimensions (M, N ) is split into tiles, resulting in MT N tiles of computation. However, as S such tiles can 2 N be processed at once simultaneously, it would result in M ST 2 batches of computation. This is depicted in Fig.8(b), where the execution time is profiled for varying tile sizes across benchmark datasets for a fixed value of T = 32.
12
Fig. 6. Total Execution Time across Benchmark Datasets profiled across all Platforms
C. Power Analysis The power dissipated by the accelerator is different for varied values of parallelism index, S, and tile size, T . For the design space exploration on power consumption, the accelerator configuration synthesized is floorplanned. The floorplanning is such that the LHS components are assigned to one block, and each of the RHS components are assigned their separate blocks. This floorplanning ensures that the design is neatly partitioned. Fig.9(a) illustrates the power breakdown of MANOJAVAM as the systolic array tile size T is increased from 4 to 20, while keeping the degree of parallelism S fixed at 4. A clear upward trend is observed across all power components—clock, signal, logic, BRAM, and DSP—with signal power exhibiting the steepest growth, followed by logic and DSP power. This trend is architecturally driven by the quadratic increase in the number of multiply-accumulate (MAC) units per systolic array, which grows as T 2 . As the tile size increases, the systolic array feeds in more routing interconnects into it, leading to higher switching activity and capacitive loads, especially along longer signal routes. This results in proportional increase in power dissipation. Logic power increases due to the expansion of operator feeding logic and controller FSMs. Concurrently, DSP power scales due to instantiation and operation of more MAC units. Interestingly, BRAM power consumption shows a marginal increase as the on chip memory footprint per array remains modest, with no additional banks instantiated as T scales. Overall, this analysis highlights the energy implications of scaling systolic array granularity, with larger tiles offering computational benefits at the cost of power and thermal overhead.
Fig.9(b) presents the impact of scaling parallelism index S on power consumption, keeping tile size constant at T = 4. As S increases from 8 to 24, a steady rise is observed across all power components, especially signal, DSP and logic power. This trend is expected, as an increment in S would instantiate a new T xT systolic array, T xT matrix accumulator, a private RHS cache unit and other interfacial components. DSP power increases due to the linear addition of MAC units. Signal power grows substantially due to the overhead of simultaneously feeding operands into more parallel datapaths, requiring wider operand buses and more frequent switching. Logic power rises as well, reflecting the growth in per-array control logic and operand routing infrastructure. Clock power shows a modest linear increase, which is consistent with the growing number of flip-flops and MAC units receiving the global clock. The BRAM power remains almost constant across the range of S due to fixed tile size. The irregularities observed at S = 12 and S = 14 is due to increased routing congestion and synthesis tool optimizations respectively. Thus, the analysis in Fig.9(b) reveals that the power dissipation increases with increase in S. D. Resource Consumption Fig.10(a) presents the FPGA resource usage profile of the accelerator, as the tile size T is increased from 4 to 20 at S = 4. All key resources - LUTs, Flip-Flops (FFs), BRAMs, and DSPs, show a monotonic increase with growing tile size. This growth is attributed to the quadratic scaling of computational units: each T xT systolic array contains T 2 MACs, leading to proportionate increase in DSP units. This scaling also drives the increase in LUTs and FFs, required
13
Fig. 7. Energy consumption across benchmarks (Log Scale). MANOJAVAM achieves up to 106 × higher efficiency than the NVIDIA A6000 by eliminating GPU driver overhead and utilizing a deterministic, low-power systolic pipeline.
(a)
(b)
Fig. 9. Design Space Exploration of Architectural Latency: (a) Impact of Tile Size (T ) on total execution time at a fixed parallelism of S = 4; (b) Impact of Parallelism (S) on total execution time at a fixed tile size of T = 4.
for operand buffering, pipeline registers and routing logic. The dip at T = 6 is associated with the synthesis tool optimization logic. BRAM utilization increases steadily, as larger tiles require deeper memory for operand storage in the caches. Fig.10(b) reflects the resource scaling of the accelerator with growing parallelism at T = 4. The resources grow linearly due to incremental instantiation of 4x4 systolic arrays, caches, accumulators and required control logic with an increase in parallelism. The DSP usage scales predictability, as 16 MAC units get instantiated with incremental rise in S. LUT and
Flip-Flop usage exhibit near identical scaling trends. BRAMs also increase gradually, with each array associated to private caches with operand B. Unlike tile size scaling, this form of spatial parallelism introduces isolated compute islands, resulting in more modular resource growth with minimal synthesis irregularities. IX. C ONCLUSION AND F UTURE S COPE OF W ORK This paper presented MANOJAVAM, a domain-specific architecture tailored to accelerate the computationally intensive
14
(a)
(b)
Fig. 10. Power Dissipation Scaling Analysis: (a) Sensitivity of power consumption to Tile Size (T ) at a fixed parallelism of S = 4; (b) Impact of Parallelism (S) on power dissipation at a fixed tile size of T = 4
(a)
(b)
Fig. 11. Hardware Resource Utilization Scaling: (a) Analysis of FPGA resource requirements (LUTs/DSPs) relative to Tile Size (T ) at S = 4; (b) Impact of Parallelism (S) on resource consumption at a fixed tile size of T = 4.
tasks of matrix multiplication and singular value decomposition (SVD) in Principal Component Analysis (PCA). The architecture introduces a novel block-streaming approach combined with a systolic array-based matrix multiplication engine and a highly pipelined Jacobian unit for executing the Jacobi algorithm efficiently. The dual-phase memory hierarchy, featuring mode-aware cache policies, further enhances data throughput and minimizes latency. MANOJAVAM has been successfully realized on both FPGA and ASIC platforms, showcasing significant speedups and power savings compared to leading CPU and GPU implementations. There are several promising directions to extending this work. MANOJAVAM can be integrated with downstream machine learning accelerators to realize a complete on-chip analytical engine for edge-AI platforms. The current architecture can also be extended to perform adaptive covariance matrix
computation to support incremental PCA. Finally, MANOJAVAM can be deployed on real life applications, such as autonomous vehicles, bio-informatics and surveillance systems to yield insights on the edge. R EFERENCES [1] R. H. Dennard, F. H. Gaensslen, H.-N. Yu, V. L. Rideout, E. Bassous, and A. R. LeBlanc, “Design of ion-implanted mosfet’s with very small physical dimensions,” IEEE Journal of solid-state circuits, vol. 9, no. 5, pp. 256–268, 1974. [2] C. A. Mack, “Fifty years of moore’s law,” IEEE Transactions on semiconductor manufacturing, vol. 24, no. 2, pp. 202–207, 2011. [3] J. Cong, M. A. Ghodrat, M. Gill, B. Grigorian, K. Gururaj, and G. Reinman, “Accelerator-rich architectures: Opportunities and progresses,” in Proceedings of the 51st annual design automation conference, 2014, pp. 1–6. [4] K. Fatahalian, J. Sugerman, and P. Hanrahan, “Understanding the efficiency of gpu algorithms for matrix-matrix multiplication,” in Proceedings of the ACM SIGGRAPH/EUROGRAPHICS Conference
15
on Graphics Hardware, ser. HWWS ’04. New York, NY, USA: Association for Computing Machinery, 2004, p. 133–137. [Online]. Available: https://doi.org/10.1145/1058129.1058148 [5] C. Silvano, D. Ielmini, F. Ferrandi, L. Fiorin, S. Curzel, L. Benini, F. Conti, A. Garofalo, C. Zambelli, E. Calore et al., “A survey on deep learning hardware accelerators for heterogeneous hpc platforms,” arXiv preprint arXiv:2306.15552, 2023. [6] H. Esmaeilzadeh, E. Blem, R. St. Amant, K. Sankaralingam, and D. Burger, “Dark silicon and the end of multicore scaling,” in Proceedings of the 38th annual international symposium on Computer architecture, 2011, pp. 365–376. [7] T. Mohaidat and K. Khalil, “A survey on neural network hardware accelerators,” IEEE Transactions on Artificial Intelligence, vol. 5, no. 8, pp. 3801–3822, 2024. [8] Q. Du and J. E. Fowler, “Hyperspectral image compression using jpeg2000 and principal component analysis,” IEEE Geoscience and Remote sensing letters, vol. 4, no. 2, pp. 201–205, 2007. [9] C. Rodarmel and J. Shan, “Principal component analysis for hyperspectral image classification,” Surveying and Land Information Science, vol. 62, no. 2, pp. 115–122, 2002. [10] A. Baobaid, M. Meribout, V. K. Tiwari, and J. P. Pena, “Hardware accelerators for real-time face recognition: A survey,” IEEE Access, vol. 10, pp. 83 723–83 739, 2022. [11] R. Pang, B. J. Lansdell, and A. L. Fairhall, “Dimensionality reduction in neuroscience,” Current Biology, vol. 26, no. 14, pp. R656–R660, 2016. [12] G. Abraham and M. Inouye, “Fast principal component analysis of largescale genome-wide data,” PloS one, vol. 9, no. 4, p. e93766, 2014. [13] K. Tsuyuzaki, H. Sato, K. Sato, and I. Nikaido, “Benchmarking principal component analysis for large-scale single-cell rna-sequencing,” Genome biology, vol. 21, no. 1, p. 9, 2020. [14] M. Greenacre, P. J. Groenen, T. Hastie, A. I. d’Enza, A. Markos, and E. Tuzhilina, “Principal component analysis,” Nature Reviews Methods Primers, vol. 2, no. 1, p. 100, 2022. [15] J. Shlens, “A tutorial on principal component analysis,” arXiv preprint arXiv:1404.1100, 2014. [16] R. Bro and A. K. Smilde, “Principal component analysis,” Analytical methods, vol. 6, no. 9, pp. 2812–2831, 2014. [17] U. A. Korat and A. Alimohammad, “A reconfigurable hardware architecture for principal component analysis,” Circuits, Systems, and Signal Processing, vol. 38, pp. 2097–2113, 2019. [18] D. Fernandez, C. Gonzalez, D. Mozos, and S. Lopez, “Fpga implementation of the principal component analysis algorithm for dimensionality reduction of hyperspectral images,” Journal of Real-Time Image Processing, vol. 16, pp. 1395–1406, 2019. [19] M. A. Mansoori and M. R. Casu, “High level design of a flexible pca hardware accelerator using a new block-streaming method,” Electronics, vol. 9, no. 3, p. 449, 2020. [20] A. Das, D. Nguyen, J. Zambreno, G. Memik, and A. Choudhary, “An fpga-based network intrusion detection architecture,” IEEE Transactions on Information Forensics and Security, vol. 3, no. 1, pp. 118–132, 2008. [21] S. N. Shahrouzi and D. G. Perera, “Dynamic partial reconfigurable hardware architecture for principal component analysis on mobile and embedded devices,” EURASIP Journal on Embedded Systems, vol. 2017, pp. 1–18, 2017. [22] T. Wu, W. Zhao, H. Guo, H. H. Lim, and Z. Yang, “A streaming pca vlsi chip for neural data compression,” IEEE transactions on biomedical circuits and systems, vol. 11, no. 6, pp. 1290–1302, 2017. [23] W. Lemaire, E. R. Koleibi, T. Omrani, M. Benhouria, K. Koua, C. Quesnel, L.-P. Gauthier, J. Ménard, K. Gagnon, S. Roy et al., “Preliminary results from a 49-channel neural recording asic with embedded spike compression in 28 nm cmos,” in 2022 20th IEEE Interregional NEWCAS Conference (NEWCAS). IEEE, 2022, pp. 285–289. [24] C. S. S. Prasanna, N. Sudha, and V. Kamakoti, “A principal component neural network-based face recognition system and asic implementation,” in 18th International Conference on VLSI Design held jointly with 4th International Conference on Embedded Systems Design. IEEE, 2005, pp. 795–798. [25] A. Elrharras, S. El Moukhlis, R. Saadane, M. Wahbi, and A. Hamdoun, “Fpga-based fully parallel pca-ann for spectrum sensing,” Computer and Information Science, vol. 8, no. 1, p. 108, 2015. [26] T. Tongyoo and Y. Ariyakul, “Fpga-based odor classification system using principal component analysis,” in 2018 International Conference on Engineering, Applied Sciences, and Technology (ICEAST). IEEE, 2018, pp. 1–4. [27] A. A. S. Ali, A. Amira, F. Bensaali, and M. Benammar, “Hardware pca for gas identification systems using high level synthesis on the zynq soc,”
in 2013 IEEE 20th International Conference on Electronics, Circuits, and Systems (ICECS). IEEE, 2013, pp. 707–710. [28] Y. Ma and D. Wang, “Accelerating svd computation on fpgas for dsp systems,” in 2016 IEEE 13th International Conference on Signal Processing (ICSP). IEEE, 2016, pp. 487–490. [29] Y.-L. Chen, C.-Z. Zhan, T.-J. Jheng, and A.-Y. Wu, “Reconfigurable adaptive singular value decomposition engine design for high-throughput mimo-ofdm systems,” IEEE transactions on very large scale integration (VLSI) systems, vol. 21, no. 4, pp. 747–760, 2012. [30] X. Wang and J. Zambreno, “An fpga implementation of the hestenesjacobi algorithm for singular value decomposition,” in 2014 IEEE International Parallel & Distributed Processing Symposium Workshops. IEEE, 2014, pp. 220–227. [31] S. Zhang, X. Tian, C. Xiong, J. Tian, and D. Ming, “Fast implementation for the singular value and eigenvalue decomposition based on fpga,” Chinese Journal of Electronics, vol. 26, no. 1, pp. 132–136, 2017. [32] M. V. Athi, S. R. Zekavat, and A. A. Struthers, “Real-time signal processing of massive sensor arrays via a parallel fast converging svd algorithm: Latency, throughput, and resource analysis,” IEEE Sensors Journal, vol. 16, no. 8, pp. 2519–2526, 2016. [33] M. U. Torun, O. Yilmaz, and A. N. Akansu, “Fpga, gpu, and cpu implementations of jacobi algorithm for eigenanalysis,” Journal of Parallel and Distributed Computing, vol. 96, pp. 172–180, 2016. [34] R. P. Brent and F. T. Luk, “The solution of singular-value and symmetric eigenvalue problems on multiprocessor arrays,” SIAM Journal on Scientific and Statistical Computing, vol. 6, no. 1, pp. 69–84, 1985. [Online]. Available: https://doi.org/10.1137/0906007 [35] R. Andraka, “A survey of cordic algorithms for fpga based computers,” in Proceedings of the 1998 ACM/SIGDA Sixth International Symposium on Field Programmable Gate Arrays, ser. FPGA ’98. New York, NY, USA: Association for Computing Machinery, 1998, p. 191–200. [Online]. Available: https://doi.org/10.1145/275107.275139 [36] G. Golub and W. Kahan, “Calculating the singular values and pseudoinverse of a matrix,” Journal of the Society for Industrial and Applied Mathematics, Series B: Numerical Analysis, vol. 2, no. 2, pp. 205–224, 1965. [37] R. Bhattacharjya, A. Sarkar, B. Maity, and N. Dutt, “Music-lite: Efficient music using approximate computing: An ofdm radar case study,” IEEE Embedded Systems Letters, vol. 16, no. 4, pp. 329–332, 2024. [38] J. E. Volder, “The cordic trigonometric computing technique,” IRE Transactions on electronic computers, no. 3, pp. 330–334, 1959. [39] P. K. Meher, J. Valls, T.-B. Juang, K. Sridharan, and K. Maharatna, “50 years of cordic: Algorithms, architectures, and applications,” IEEE Transactions on Circuits and Systems I: Regular Papers, vol. 56, no. 9, pp. 1893–1907, 2009. [40] J. Zhuang, J. Lau, H. Ye, Z. Yang, S. Ji, J. Lo, K. Denolf, S. Neuendorffer, A. Jones, J. Hu et al., “Charm 2.0: Composing heterogeneous accelerators for deep learning on versal acap architecture,” ACM Transactions on Reconfigurable Technology and Systems, vol. 17, no. 3, pp. 1–31, 2024. [41] C. S. Mummidi, V. C. Ferreira, S. Srinivasan, and S. Kundu, “Highly efficient self-checking matrix multiplication on tiled amx accelerators,” ACM Transactions on Architecture and Code Optimization, vol. 21, no. 2, pp. 1–22, 2024. [42] N. P. Jouppi, C. Young, N. Patil, D. Patterson, G. Agrawal, R. Bajwa, S. Bates, S. Bhatia, N. Boden, A. Borchers et al., “In-datacenter performance analysis of a tensor processing unit,” in Proceedings of the 44th annual international symposium on computer architecture, 2017, pp. 1–12. [43] A. S. Prasad, M. Scherer, F. Conti, D. Rossi, A. Di Mauro, M. Eggimann, J. T. Gómez, Z. Li, S. S. Sarwar, Z. Wang et al., “Siracusa: A 16 nm heterogenous risc-v soc for extended reality with at-mram neural engine,” IEEE Journal of Solid-State Circuits, 2024. [44] Y. Zhang, X. Zhang, P. Xu, Y. Zhao, C. Hao, D. Chen, and Y. Lin, “Autoai2c: An automated hardware generator for dnn acceleration on both fpga and asic,” IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, 2024. [45] H. Chen, Y. Ni, A. Zakeri, Z. Zou, S. Yun, F. Wen, B. Khaleghi, N. Srinivasa, H. Latapie, and M. Imani, “Hdreason: Algorithm-hardware codesign for hyperdimensional knowledge graph reasoning,” arXiv preprint arXiv:2403.05763, 2024. [46] W. Ren, S. Koteshwara, M. Ye, H. Franke, and D. Chen, “S2tar: Shared secure trusted accelerators with reconfiguration for machine learning in the cloud,” in 2024 IEEE 17th International Conference on Cloud Computing (CLOUD). IEEE, 2024, pp. 267–278. [47] Y. Xue, Y. Liu, L. Nai, and J. Huang, “Hardware-assisted virtualization of neural processing units for cloud platforms,” in 2024 57th IEEE/ACM
16
International Symposium on Microarchitecture (MICRO). IEEE, 2024, pp. 1–16. [48] S. Mondal, S. D. Manasi, K. Kunal, Z. Zeng, S. S. Sapatnekar et al., “A unified engine for accelerating gnn weighting/aggregation operations, with efficient load balancing and graph-specific caching,” IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, vol. 42, no. 12, pp. 4844–4857, 2022. [49] B. Zhang and V. Prasanna, “Dynasparse: Accelerating gnn inference through dynamic sparsity exploitation,” in 2023 IEEE International Parallel and Distributed Processing Symposium (IPDPS). IEEE, 2023, pp. 233–244. [50] K. Marino, P. Zhang, and V. K. Prasanna, “Me-vit: A single-load memory-efficient fpga accelerator for vision transformers,” in 2023 IEEE 30th International Conference on High Performance Computing, Data, and Analytics (HiPC). IEEE, 2023, pp. 213–223. [51] K.-Y. Chen, C.-S. Yang, Y.-H. Sun, C.-W. Tseng, M. Fayazi, X. He, S. Feng, Y. Yue, T. Mudge, R. Dreslinski et al., “Dap: A 507-gmacs/j 256-core domain adaptive processor for wireless communication and linear algebra kernels in 12-nm finfet,” IEEE Journal of Solid-State Circuits, 2024. [52] Q. Zhang, Z. Fan, H. An, Z. Wang, Z. Li, G. Wang, P. Abillama, H.-S. Kim, D. Blaauw, and D. Sylvester, “Robovisio: A micro-robot vision domain-specific soc for autonomous navigation enabling fully-on-chip intelligence via 2-mb emram,” IEEE Journal of Solid-State Circuits, 2024. [53] J.-F. Zhang, C.-H. Lu, and Z. Zhang, “Tetrix: Flexible architecture and optimal mapping for tensorized neural network processing,” IEEE Transactions on Computers, 2024. [54] F. H. McMahon, “The livermore fortran kernels: A computer test of the numerical performance range,” Lawrence Livermore National Lab., CA (USA), Tech. Rep., 1986. [55] N. P. Jouppi, “Cache write policies and performance,” ACM SIGARCH Computer Architecture News, vol. 21, no. 2, pp. 191–201, 1993. [56] J. Dongarra, C. Moler, and R. James, “Bunch, and gw stewart,” LINPACK users’ guide. SIAM, 1979. [57] I. Bravo, M. Mazo, J. L. Lázaro, A. Gardel, P. Jiménez, and D. Pizarro, “An intelligent architecture based on field programmable gate arrays designed to detect moving objects by using principal component analysis,” Sensors, vol. 10, no. 10, pp. 9232–9251, 2010. [58] S. Kasap and S. Redif, “Novel field-programmable gate array architecture for computing the eigenvalue decomposition of para-hermitian polynomial matrices,” IEEE Transactions on Very Large Scale Integration (VLSI) Systems, vol. 22, no. 3, pp. 522–536, 2013. [59] E. Alpaydin and C. Kaynak, “Optical Recognition of Handwritten Digits,” UCI Machine Learning Repository, 1998, DOI: https://doi.org/10.24432/C50P49. [60] Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner, “Gradient-based learning applied to document recognition,” Proceedings of the IEEE, vol. 86, no. 11, pp. 2278–2324, 1998. [61] A. Krizhevsky, G. Hinton et al., “Learning multiple layers of features from tiny images,” 2009. [62] “Olivetti faces dataset,” AT&T Laboratories Cambridge (1992–1994), available via scikit-learn, 1994, https://scikit-learn.org/0.19/datasets/ olivetti faces.html. [63] W. Wolberg, O. Mangasarian, N. Street, and W. Street, “Breast Cancer Wisconsin (Diagnostic),” UCI Machine Learning Repository, 1993, DOI: https://doi.org/10.24432/C5DW2B. [64] T. Mitchell, “Twenty Newsgroups,” UCI Machine Learning Repository, 1997, DOI: https://doi.org/10.24432/C5C323. [65] NVIDIA Corporation, “NVIDIA Ampere GA102 GPU Architecture Whitepaper,” NVIDIA Corporation, Tech. Rep., 2020. [Online]. Available: https://www.nvidia.com/content/PDF/ nvidia-ampere-ga-102-gpu-architecture-whitepaper-v2.pdf
Srivaths Ramasubramanian (Student Member, IEEE) received the B.E. degree in Electronics and Communication Engineering from Rashtreeya Vidyalaya College of Engineering (RVCE), Bengaluru, India, in 2025. He is currently pursuing the Ph.D. degree in Computer Engineering at the Donald Bren School of Information and Computer Sciences, University of California, Irvine. His research interests include hardware acceleration for machine learning and communication systems, computer architecture and FPGA-based system design.
Anjali Devarajan (Student Member, IEEE) received the B.E. degree in Electronics and Communication Engineering from Rashtreeya Vidyalaya College of Engineering (RVCE), Bengaluru, India, in 2025. She is currently an engineer at Infineon Technologies. Her research interests are embedded systems, VLSI design, and hardware acceleration for machine learning. Her experience includes internships in ADAS, control systems and low-power design. Her work has been published in IEEE and Springer conferences.
Kousthub P Kaivar (Student Member, IEEE) received the B.E. degree in Electronics and Communication Engineering from Rashtreeya Vidyalaya College of Engineering (RVCE), Bengaluru, India, in 2025. He is currently a Hardware Engineer in the Battery Management Unit team at Enphase Energy. His interests include VLSI, Computer Architecture, FPGA based system design, communication systems and machine learning. In addition, he has work published in communication systems involving antennas and radars in IEEE conferences.
Vibha Shrestta (Student Member, IEEE) received the B.E. degree in Electronics and Communication Engineering from Rashtreeya Vidyalaya College of Engineering (RVCE), Bengaluru, India, in 2025. She is currently working as a Design engineer at ABB. Previously, she has also worked in the research and development wing at DRDO, solving critical challenges in embedded system design and flight avionics. Her interests include VLSI, FPGA based system design, analog electronics and power electronics.
Shashank D (Student Member, IEEE) received the B.E. degree in Electronics and Communication Engineering from Rashtreeya Vidyalaya College of Engineering (RVCE), Bengaluru, India, in 2025. He is currently working as an Associate Sales Development Representative at Astar Data LLP . His interests include VLSI, FPGA based system design, communication systems and Machine Learning.
17
Dr. Sowmyarani C.N (Member, IEEE) received the Ph.D. degree in data privacy. She is currently a Professor with the Department of Computer Science and Engineering, R. V. College of Engineering. She has 16 years of teaching experience. She has published more than 60 research publications. She has delivered many talks and training sessions in the field of data privacy, networks, and security. She has been working on research projects funded by VGST, Government of Karnataka, India. Her research interests include data engineering, data privacy, computer networks, cyber security and education technology.
Dr. Govinda Raju M received the B.E. degree in Electronics and Communication Engineering from Visveswaraya Technological University(VTU), Karnataka, India, in 2004, and the M.Tech. and Ph.D. degrees in Electronics from VTU in 2007 and 2013, respectively. He is currently an associate professor in the Department of Electronics and Communication Engineering at RV College of Engineering, Bengaluru, India. His research interests include embedded system design, computer architecture, real-time systems, and embedded automotive systems. He has published over 50 research papers in reputed journals and conferences. He has guided several postgraduate and undergraduate student projects funded by national bodies.
Dr. K.S Geetha (Senior Member, IEEE) received the B.E. degree in Electronics and Communication Engineering in 1991 and the M.Tech. degree in Computer Applications to Industrial Drives in 1998 from the National Institute of Engineering, Mysore, India, and the Ph.D. degree from Visvesvaraya Technological University, India, in 2012. She is a professor in the department of Electronics & Communication Engineering and currently the Vice Principal of Rashtreeya Vidyalaya College of Engineering (RVCE), Bengaluru, India. Her research interests include image processing, signal processing, VLSI design, and flexible electronics.