Reproducible Ozaki-Based Matrix Multiplication on Photonic Matrix Processors With Algorithm-Based Fault Tolerance

2.1 Photonic Matrix Processors

Photonic matrix multiplication processors can be classified into three architectures: free-space optical systems, integrated photonic circuits using Mach-Zehnder interferometer (MZI) meshes, and wavelength-division multiplexing (WDM)-based approaches [8]. Among these, MZI-based integrated photonic processors have attracted significant attention due to their scalability and compatibility with silicon photonics platforms.

Optical AI inference accelerators are attracting attention [1, 2]. Recently, accelerators with photonic tensor cores, RISC-V cores, and embedded OSs have appeared, and they can perform multiplications of $$$ 64\times 64 $$$ matrices and $$$ 64\times 1 $$$ matrices, and $$$ 128\times 128 $$$ matrices and $$$ 128\times 1 $$$ matrices at once for numerical values with about 8-bit accuracy [2]. It is fabricated on GlobalFoundries' 45 nm silicon photonics platform, while advanced CMOS nodes (e.g., 12 nm) are used externally for high-performance electronic control/processing. In this way, research on photonic accelerators is moving to a level where contributions can be made not only from optical science and physics, but also from the fields of computer science.

Despite their high throughput and energy efficiency, photonic matrix processors inherently suffer from limited numerical precision due to analog noise sources such as shot noise, thermal noise, and device non-idealities. This limitation restricts their applicability to low-precision workloads such as AI inference, and creates a fundamental challenge for extending their use to high-precision scientific computing.

Another concern is the matrix-reconfiguration overhead of photonic matrix processors. While reconfiguration latency can be a concern in current implementations, it is expected to diminish with advances in device technologies such as fast tunable photonic elements, integrated control circuits, and improved calibration techniques. Moreover, emerging use cases, particularly in AI workloads with repeated or slowly changing weight matrices, naturally amortize reconfiguration costs over many computations. Therefore, although matrix reconfiguration remains an important system-level consideration, it is orthogonal to the core contribution of this study and is expected to become a less dominant factor in future photonic computing platforms.

2.2 Precision Limitations in Analog Computing

The precision limitation of analog computing systems has been widely studied in the context of compute-in-memory (CIM) and photonic computing. Commercial CIM products have recently emerged. Mythic has developed analog in-memory computing processors that perform matrix operations within memory arrays for low-power AI inference [9]. In CIM architectures, dominant error sources include device variability, nonlinearity, and IR drop, which often manifest as deterministic or structured errors that can be mitigated through calibration and device-aware mapping [10].

In contrast, photonic matrix processors are dominated by stochastic noise sources that can be modeled as additive perturbations, such as additive white Gaussian noise (AWGN). Prior work has shown that such noise fundamentally limits achievable floating-point precision, particularly in accumulation-heavy operations such as matrix multiplication [8]. As a result, there exists a practical precision ceiling beyond which increasing arithmetic precision does not improve effective computational accuracy.

2.3 High-Precision Computation With Low-Precision Units

To address precision limitations, several techniques have been proposed to enable high-precision computation using low-precision hardware. Mixed-precision methods, such as those used in GPU tensor cores, combine low-precision arithmetic with iterative refinement to recover high-precision results [11].

The Ozaki scheme is a representative approach that decomposes high-precision matrix multiplication into multiple low-precision operations [4, 5]. The original description of the scheme commonly represents each slice in a low-precision floating-point format. This study instead targets the integer-based implementation of the Ozaki scheme in [5]. It divides the shared mantissa space of each row of $$$ A $$$ and each column of $$$ B $$$ into successive bit ranges and stores the extracted mantissa slices as $$$ \mathrm{INT}Q $$$ matrices.

This representation is related to block floating point because all integer values in one row slice of $$$ A $$$, or one column slice of $$$ B $$$, share a power-of-two exponent. The shared exponents are retained as metadata by the xPU, whereas only the integer mantissa slices are sent to the photonic matrix-multiplication unit. Thus, block scaling here refers to this shared-exponent representation of each integer Ozaki slice. Microscaling is another possible shared-scale representation [12], but it is not assumed in this study.

However, existing applications of the Ozaki scheme do not explicitly consider stochastic analog noise on matrix-multiplication hardware units.

2.4 Algorithm-Based Fault Tolerance

In digital system communications, as a method of dealing with soft errors (temporary bit corruption), the sender adds an ECC (Error Checking and Correcting) code to the original data before transferring it, allowing the receiver to detect and correct soft errors within a certain number of errors.

Similarly, in ABFT, an array for error detection and correction is added to the original array in advance in matrix multiplication, making it possible to detect and correct soft errors within a certain number of errors from the calculation results [6, 13, 14]. ABFT is characterized by its ability to provide reliability for the entire hardware process, including arithmetic units and memory access, because it operates at the application level. In floating-point ABFT, distinguishing soft errors from rounding errors requires carefully designed detection thresholds. Prior work has explored adaptive and distributed threshold estimation, such as A-ABFT and V-ABFT, to control false positives and false negatives [15]. In contrast, our target is integer ABFT.

Most existing ABFT techniques are designed for discrete error models, such as bit flips in digital systems. Their applicability to analog computing environments, where errors are continuous-valued and stochastic, remains largely unexplored. In particular, it is unclear how ABFT interacts with floating-point representations under analog noise.

2.5 Gap and Contribution

Despite extensive research on precision recovery and fault tolerance, there has been limited work on integrating high-precision decomposition techniques with algorithm-level error correction under stochastic analog noise. Existing approaches either focus on deterministic error models or rely on calibration-based mitigation techniques that are not directly applicable to photonic systems.

In this paper, we address this gap by combining the Ozaki scheme with the ABFT under an AWGN-based noise model. Our approach provides a unified framework for enabling high-precision floating-point matrix multiplication on low-precision photonic processors, explicitly considering the impact of stochastic analog noise on floating-point representations.

3.1 Computation Model and Decomposition

We consider a photonic tensor core that performs matrix-vector multiplication of size $$$ N\times N $$$ with an $$$ N $$$-dimensional vector in a single operation. The overall future 2.5D integrated system is illustrated in Figure 1, which abstracts a photonic AI acceleration processor, which co-packaged a silicon photonics die and CMOS drivers, such as [2].

Our goal is to compute matrix-matrix multiplication $$$ \mathbf{C}=\mathbf{AB} $$$, where $$$ \mathbf{A},\mathbf{B}\in {\mathbb{R}}^{M\times M} $$$, using such tensor cores.

This decomposition indicates that computation is dominated by repeated matrix-vector operations.

3.2 Numerical Precision Recovery via Ozaki Decomposition

Photonic tensor-cores operate at low precision due to analog noise and limited bit representation. To enable high-precision computation, we employ the integer-based Ozaki scheme in [5], which decomposes high-precision matrix multiplication into multiple low-precision operations.

Each submatrix multiplication $$$ {\mathbf{A}}_i{\mathbf{B}}_j $$$ is realized as a sequence of matrix-vector multiplications using the photonic tensor core, while accumulation is performed on a digital processor (xPU). This decomposition allows arbitrary precision by increasing $$$ s $$$ and $$$ t $$$, but also increases the number of partial computations, making the system more vulnerable to noise-induced errors.

3.3 ABFT for Error Detection and Correction

To ensure reliability under stochastic analog noise, we introduce a simple ABFT [6] at the level of each decomposed computation.

Due to the linearity of matrix multiplication, checksum consistency is preserved. Errors can be detected by comparing the checksum values with recomputed sums. Furthermore, the error location can be identified as the intersection of inconsistent row and column checksums, enabling single-error correction. This holds regardless of whether the error occurs in a data element or in a checksum element, due to the orthogonal redundancy provided by row and column checksums.

We illustrate the ABFT using a simple $$$ 2\times 2 $$$ matrix multiplication.

Thus, the erroneous element is identified as the intersection of the inconsistent row and column.

This example demonstrates that the ABFT can both detect and correct a single error within a protection unit. The same principle holds for a general $$$ N\times N $$$ protection unit. Assume that checksum consistency holds in the error-free case, and that at most one element in a protection unit is corrupted. Then, the erroneous element can be uniquely identified as the intersection of the inconsistent row and the inconsistent column, and its correct value can be recovered from the corresponding checksum relation. If two or more elements are corrupted within the same protection unit, the inconsistent rows and columns may no longer determine a unique intersection, and correct recovery is not guaranteed. This example and the above proposition demonstrate that the ABFT can both detect and correct a single error within a protection unit. It also differs from randomized detection methods, that is, Freivalds' algorithm, in that the recovery is deterministic within a protection unit and does not require re-execution after error detection.

3.4 Procedure

This procedure ensures that high-precision results are obtained from low-precision photonic computations and xPU's preprocessing and postprocessing. The original input matrices and the final reconstructed result are represented in floating-point format, that is, FP64 in this study. The xPU determines the shared scale vectors and extracts the integer mantissa slices by bit masking. Only these $$$ \mathrm{INT}Q $$$ slices are processed by the photonic tensor core. Checksum verification and error correction are performed using exact equality comparison on the integer matrices and their integer outputs. After error correction, the xPU applies the power-of-two row and column scales and accumulates the partial products in FP64.

The xPU stores the FP64 input matrices, the integer slices, the shared scale vectors, and the FP64 accumulator. It performs scale selection, mantissa-bit extraction, checksum generation and verification, single-error correction, power-of-two rescaling, and accumulation. The scale metadata is small compared with the matrix data, whereas the main xPU arithmetic costs are slice generation, checksum processing, and floating-point accumulation of the partial products.

4.1 Evaluation Objectives and Assumptions

The goal of this evaluation is to understand how stochastic analog noise affects reliability, and how the proposed method mitigates these effects.

These analyses collectively reveal the design tradeoffs between precision, redundancy, and noise tolerance.

The parameter values used in this evaluation are shown in Table 1. Each value in the target matrix multiplication is a double-precision floating-point number (FP64). Furthermore, the photonic tensor core supports $$$ Q $$$-bit number representation, that is, $$$ \mathrm{INT}Q $$$, for input values, and $$$ 3Q $$$-bit representation for output in order to prevent overflow in the multiplication and sum results. The output of the photonic MAC requires a wider bit width than the input precision due to accumulation. We assume that the Ozaki scheme works with the fast mode. The fast mode of the Ozaki scheme reduces computational complexity from $$$ O(st) $$$ to approximately $$$ \frac{s\left(s+1\right)}{2} $$$ by truncating low-order components and cross terms, trading exact reconstruction for efficiency at the cost of structured truncation errors. The value of $$$ s $$$ is picked up from the results in [16].

In the GMM-based noise model, we set the burst probability $$$ {\pi}_B=0.01 $$$ and the variance scaling factor $$$ \alpha =10 $$$, where the high-noise state is defined as $$$ {\sigma}_B=\alpha {\sigma}_G $$$. These values are chosen as representative parameters to capture burst-like error behavior observed in analog and photonic computing systems, where rare but significantly larger noise events occur intermittently. The chosen $$$ {\pi}_B $$$ reflects a low-probability burst regime, while $$$ \alpha $$$ represents a moderate-to-severe noise amplification during burst events. Rather than targeting a specific hardware instance, these parameters are intended to provide a qualitative understanding of how correlated noise affects the reliability of the proposed method.

4.2 Computational Cost Analysis

When performing FP64 matrix operations (multiply-and-accumulate calculations) using the proposed method, we counted the number of calls to the photonic tensor core, taking into account this ABFT, in addition to the number of calls to the integer matrix operations in [5].

This produces staircase-like curves, especially around tiling boundaries.

The number of calls to the photonic tensor core is approximated in Figures 2 and 3. Note that in digital circuits, the operations required for the Ozaki scheme can be simply realized with shift and addition operations when the original matrix is FP64.

In Figures 2 and 3, the horizontal axis represents the FP64 input matrix size ($$$ M $$$), and the vertical axis represents the number of invocations of the photonic tensor core. ABFT performs redundant calculations regardless of whether errors exist or not. Therefore, the number of invocations does not depend on the analog error levels.

Furthermore, as expected, we can see that the number of invocations of the photonic tensor core is determined by the matrix size that can be executed at one time ($$$ N $$$).

4.3 ULP Error Analysis

We analyze the impact of stochastic analog noise on integer quantized outputs produced by photonic tensor cores. We consider two noise models: AWGN and a two-state Gaussian mixture (burst-like) for representing independent and identically distributed and burst-like behavior via mixtures of low- and high-variance Gaussian components, respectively.

4.3.3 Analytical Results

In Figure 4, the curves for INT4, INT8, and INT16 for AWGN exhibit identical shapes but are horizontally shifted. This behavior arises because the error probability is governed by the normalized ratio. Therefore, increasing the bit width effectively shifts the tolerance to smaller noise magnitudes by an exponential factor. The analysis highlights a fundamental precision and noise tradeoff in analog photonic computation. While increasing the bit width improves numerical resolution, it simultaneously reduces robustness to stochastic noise. Importantly, although the ULP is constant ($$$ =1 $$$) in the integer code domain, the effective ULP in the normalized real-value domain decreases with bit width, which ultimately determines error behavior. As expected, this distinction explains why integer representations can exhibit strong bit-width-dependent sensitivity under stochastic analog noise, despite having a constant ULP in the code domain.

Another concern is the comparison of AWGN cases with the GMM cases. Unlike AWGN, the GMM introduces temporally correlated high-noise states, increasing the probability of multiple simultaneous errors within a protection unit, which directly undermines the single-error correction assumption of ABFT. Under AWGN, all samples are drawn from a single moderate-variance distribution, leading to consistently higher error probabilities. In contrast, GMM spends most of the time in a low-variance state, reducing the average error rate despite occasional high-noise events.

Figure 5 plots the probability of more than one error for INT4, INT8, and INT16. The $$$ X $$$-axis is the noise standard deviation in the normalized domain $$$ \left[-1,1\right] $$$. Since they can recover at most one error correction, the smaller is better in $$$ Y $$$-axis for the proposed method. This result highlights a fundamental limitation of ABFT under burst-like noise: as the protection unit size increases, the probability of uncorrectable multi-error events grows rapidly, limiting the scalability of the proposed approach.

4.4 Failure Rate

Even with precision decomposition and algorithm-level error correction, residual errors may remain. In this section, we analytically characterize the failure probability of the proposed method under stochastic analog noise.

4.4.4 Analytical Results

Assume that photonic tensor core provides 8-bit precision, $$$ N=256 $$$, Figure 6 plots the failure probability without error correction. It uses the fast-mode Ozaki scheme. As expected, the failure probability increases rapidly with matrix size, and large-scale computations become unreliable even at moderate error rates. Figure 7 shows the failure probability with the ABFT. The results demonstrate a substantial reduction in failure probability across all matrix sizes, confirming the quadratic suppression effect predicted by the analytical model.

Overall, the analysis demonstrates that although the Ozaki scheme increases the number of sub-computations, the integration of ABFT enables scalable and reliable high-precision computation on noisy photonic hardware.

Figures 8 and 9 show the failure probabilities under the GMM model for the same configurations as Figures 6 and 7. The $$$ X $$$-axis represents the base noise standard deviation ($$$ \sigma $$$) of the stochastic analog noise added to the photonic MAC output before quantization. Compared to the AWGN case, the GMM results exhibit a significantly higher failure probability, particularly in the low-error regime. This is because the burst-like model increases the likelihood of multiple simultaneous errors within the same protection unit. While ABFT effectively suppresses failure probability under AWGN due to its single-error correction capability, its effectiveness degrades under GMM, where correlated errors violate the single-error assumption. As a result, even when the average noise level is low, the presence of rare but high-variance noise events leads to a sharp increase in failure probability.

These results indicate that ABFT effectiveness is strongly dependent on the error distribution, not only the average error rate. In particular, burst-like error patterns can significantly degrade reliability, even when the mean error rate remains unchanged. That is, the reduction in failure rate achieved by ABFT is significantly smaller under GMM than under AWGN.

4.5 System-Level Execution-Time Model

The execution time is influenced not only by computation and memory transfer, but also by analog-to-digital conversion (ADC) latency. In photonic tensor cores, each analog multiply-accumulate (MAC) result must be converted into a digital representation before further processing, introducing a non-negligible overhead.

In the evaluation, we assume $$$ {P}_{\mathrm{eff}}=10 $$$ TFLOPS, $$$ {B}_{\mathrm{eff}}=500 $$$ GB/s, and $$$ {t}_{\mathrm{adc}}=5 $$$ ns. The input precision is modeled with $$$ {b}_{\mathrm{in}}=1 $$$ byte, and accumulation with $$$ {b}_{\mathrm{acc}}=8 $$$ bytes. Reuse factors are set to $$$ {r}_A={r}_B=16 $$$ and $$$ {r}_C=1 $$$.

For the comparison between no ABFT and ABFT, both methods use the same Ozaki decomposition count. Therefore, the difference arises only from the additional checksum rows and columns.

Figure 10 shows the execution time for each method. Each panel separates the computation and ADC components and directly labels the total execution time above each bar. The relative ABFT overhead is also shown in each panel. It decreases from $$$ 0.587\% $$$ at $$$ M=512 $$$ to $$$ 0.073\% $$$ at $$$ M=4096 $$$. Under the evaluated parameters, computation plus ADC latency is larger than the overlapped memory-transfer time and therefore determines the total shown in the figure.

More specifically, the overhead of the ABFT decreases as $$$ M $$$ increases, since the relative impact of a single checksum row and column diminishes. These results indicate that, when the Ozaki decomposition count is fixed, the performance difference among the three methods is small. In the exact discrete tiling model, the effective payload size per tensor-core invocation becomes $$$ \left(N-1\right)\times \left(N-1\right) $$$ when ABFT is applied, due to the inclusion of checksum rows and columns within the same array. This modifies the invocation count from $$$ {M}^3/{N}^2 $$$ to approximately $$$ {M}^3/{\left(N-1\right)}^2 $$$. However, since $$$ {\left(N/\left(N-1\right)\right)}^2\approx 1 $$$ for practical values of $$$ N $$$ (e.g., 1.03 for $$$ N=64 $$$ and 1.008 for $$$ N=256 $$$), this effect introduces only a minor constant-factor overhead and does not change the overall execution-time trend.

Overall, under the model with overlap and data reuse, the system behaves as compute-bound including ADC rather than memory-bound. This suggests that conventional architectural techniques such as blocking, buffering, and pipelining are effective in mitigating memory overhead.

Notice that the matrix-reconfiguration overhead of photonic matrix processors is not explicitly considered in this study, as our focus is on the intrinsic performance potential of analog photonic MAC operations. While reconfiguration latency can be a concern in current implementations, it is expected to diminish with advances in device technologies such as fast tunable photonic elements, integrated control circuits, and improved calibration techniques. Moreover, emerging use cases, particularly in AI workloads with repeated or slowly changing weight matrices, naturally amortize reconfiguration costs over many computations. Therefore, although matrix reconfiguration remains an important system-level consideration, it is orthogonal to the core contribution of this study and is expected to become a less dominant factor in future photonic computing platforms.