






















\section*{Acknowledgement}
This work is supported by NSF IIS-2008102.

\nocite{langley00}

\section{Experiments}







A MCTensor is able to achieve improved precision compared with a PyTorch Tensor with the same data type. If higher precision implies better results under the same experiment settings, we would expect our $n$-MCTensor ($n \geq 2$) model can achieve better performance than a native PyTorch model with same precision. We demonstrate the increasing precision under the same data type for MCTensor by carrying out experiments in two settings: (1) high precision learning with low precision numbers, this include logistic regression model
and multi-layer perceptrons (MLP) in Section~\ref{subsec:mlp}; and (2) high precision hyperbolic embedding reconstruction task in Section~\ref{subsec:hyperbolic}.
For all these experiments, MCTensor models use the same initialization as the Tensor models, and use \pythoninline{MCModule} as its PyTorch's \pythoninline{nn.Module} counterpart, with MCTensors as layer weights. Our code and models are open-source and available at github \footnote{\url{https://github.com/ydtydr/MCTensor}}.

 
 

 

\subsection{High-precision Computations with Low Precision}
Low precision machine learning employs models that compute with low precision numbers (e.g. float8 and float16), which become popular in many edge-applications \cite{hubara2017quantized,zhang2019qpytorch}. However, the quantization error inherent with low precision arithmetic could affect the convergence and performance of the model \cite{wu2018training,de2018high}. As a matter of fact, most current low precision learning frameworks including \cite{courbariaux2014training,de2018high} adopt high precision numbers during gradient computations and optimizations, but pursue a better way to convert the results to low-precision numbers with less quantization errors. In comparison, with MCTensor, one can get \emph{accurate} update with purely low-precision numbers thoroughly without even touching high precision arithmetics. This is particularly helpful on devices where only low-precision arithmetics are supported.


\paragraph{Logistic Regression.}
\label{subsec:logistic}
\begin{figure}
    \centering
    \includegraphics[width=0.4\textwidth]{figs/experiments/logreg/logistic_cancer_train.png} 
    \vspace{-1em}
    \caption{\centering Training Losses for Logistic Regression Models on the Breast Cancer Dataset}
       
    \label{fig:my_label}
    \vspace{-1em}
\end{figure}

We conduct a logistic regression task on a synthetic dataset and the cancer dataset, both are datasets with binary labels. The synthetic dataset consists of $1,000$ data points, where each data point contains two features. This dataset is constructed through the \pythoninline{make\_classification} \cite{guyon2003design} function from scikit-learn package with both \pythoninline{n\_informative} and \pythoninline{n\_clusters\_per\_class} are set to 1.  The breast cancer dataset \cite{breast_cancer} consists of 569 data points and each data point contains 30 features. More details can be found in \ref{logreg}.









\paragraph{Multi-layer Perceptron.}
\label{subsec:mlp}
We also construct a Multi-layer Perceptron using MCTensor (MC-MLP) and evaluate it on classification tasks on the breast cancer dataset and a reduced MNIST, which has only 10,000 data points in total, with 1,000 images sampled randomly per class \cite{deng2012mnist}. The MC-MLP consists of three \pythoninline{MCLinear} layers with model weight in MCFloat16. After each \pythoninline{MCLinear} layer, the resulting MCTensor is transformed into a normal tensor, passed it to activation function and fed to the next layer.
\begin{figure}[t]
    \centering
    \includegraphics[width=0.4\textwidth]{figs/experiments/mlp_cancer/mlp_cancer_test.png}
    \vspace{-1.5em}
    \caption{Test Accuracy for MC-MLP on Breast Cancer Dataset. Notice that the curve for MC-MLP model with nc=2, float16 essentially overlaps with the curve for PyTorch MLP model with float32.}
    \label{fig:my_label}
\end{figure}

\begin{figure}[t]
    \centering
    \includegraphics[width=0.4\textwidth]{figs/experiments/mlp_mnist/MNIST-MLP-test.png}
    \vspace{-1.5em}
    \caption{Test Accuracy for MC-MLP on Reduced MNIST Dataset. Notice that the curve for MC-MLP model with nc=2, float16 essentially overlaps with the curve for PyTorch MLP model with float32.}
    \label{fig:mcmlp-test-mnist}
\end{figure}

We experiment on both the Breast Cancer dataset and the Reduced MNIST dataset, and as demonstrated in Table~\ref{tab:mlp-b} and Figure~\ref{fig:mcmlp-test-mnist}, in both experiments, MC-MLP outperforms PyTorch float16 models, and arrive at a lower training loss and a higher test accuracy (same as PyTorch float32 and float64 models). Notice that for both datasets, after $nc$ exceeds certain value ($nc=2$), adding extra more $nc$ would not lead to further improvement. 
More details can be found in Section \ref{MC-MLP}.

\input{tables/mlp-b}










    
    







\subsection{Hyperbolic Embedding} 
\label{subsec:hyperbolic}
Another use case that benefits from MCTensor is hyperbolic deep learning, where the non-Euclidean hyperbolic space is adopted in place of Euclidean space for various purposes. It has been shown in \cite{nickel2017poincare,chami2019hyperbolic,yu2022hyla} that hyperbolic space is better suited for processing hierarchical data (e.g. trees, acyclic graphs). However, representing hyperbolic space with ordinary floating-points (even float64) leads to unbounded representation error, particularly when points get far away from the origin, known as the ``NaN" problem \cite{yu2019numerically} that not-a-number error occurs during computations. With MCTensor, we can achieve high-precision computations and avoid the ``NaN" problem by simply using more components. 

Following \citet{yu2021representing}, we conduct hyperbolic embedding experiment on the WordNet Mammals dataset with 1181 nodes and 6541 edges. We use Poincaré upper-half space (\textbf{Halfspace}) model to embed its transitive closure for reconstruction. The Poincaré Halfspace model is the manifold $(\mathcal{U}^n, g_u)$, where $\mathcal{U}^n$ is the upper half space of the $n$-dimensional Euclidean space, with the metric tensor and distance function being:
\[
\textstyle
g_u(\mathbf{x}) = \frac{g^E}{x_n^2},~~~~ d_u(\mathbf{x},\mathbf{y}) = \text{arcosh}
\left(1+ \frac{||\mathbf{x}-\mathbf{y}||^2}{2x_ny_n} \right),
\]
where $g^E$ is the Euclidean metric tensor. For the observed edges $\mathcal{D} = \{(\mathbf{x},\mathbf{y})\}$, we learn the embeddings $\Theta$ for all nodes, subject to minimizing the following loss function
\[
\textstyle
\mathcal{L}(\Theta) = \sum_{(\mathbf{x},\mathbf{y})\in D} \log \frac{e^{-d_u(\mathbf{x},\mathbf{y})}}{\sum_{\mathbf{y'}\in \mathcal{N}(\mathbf{x})} e^{-d_u(\mathbf{x},\mathbf{y'})}},
\]
where $\mathcal{N}(\mathbf{x})$ are randomly chosen 50 negative examples in addition to the positive example $(\mathbf{x},\mathbf{y})$. We report the results in Table~\ref{tab:mammals}, where MAP is the mean averaged precision and MR is the mean rank metric.


\input{tables/Mammal}


 
\section{Introduction}

High precision computations are of interest in many areas. For example, an emerging trend in studying dynamical systems is to use Taylor methods with high-precision arithmetic \cite{bailey2015high}, and delaunay triangulation in computational geometry \cite{schirra1998robustness}. Recently, high precision computations are even desired for some deep learning tasks, e.g. hyperbolic deep learning, so as to use hyperbolic space stably in practice \cite{yu2019numerically, yu2021representing}.

As of now, most high precision algorithms can be divided in two categories: (1) the standard multiple-digit ``BigFloat" format using a sequence of digits coupled with a single exponent term; (2) the multiple-component floating-point format (MCF) using an unevaluated sum of multiple ordinary floating-point numbers (e.g. float16, float32). Some examples of the first approach include the GNU Multiple Precision (GMP) Arithmetic Library \cite{Granlund12} and Julia's BigFloat type \cite{bezanson2017julia}. The idea of multiple-component float approach (also referred as ``expansion" in some literature), dates back to priest's works \cite{priest1991algorithms}, and some example implementations include multi-component library on C++ \cite{hida2007library} and MATLAB \cite{jiang2016implementation}. The multiple-digit approach can represent compactly a much larger range of numbers, whereas the multiple-component approach still adopts the limited precision floats to achieve high precision computations with error guarantees \cite{joldecs2015arithmetic}. 
However, the multiple-component approach has an advantage in speed over the multiple-digit approach, as it uses standard floats such as float16 and float32 and hence better takes advantage of existing floating-point accelerators \cite{hida2007library,shewchuk1997adaptive}.


In the emerging deep learning area, there is a demand for high precision computations in some deep learning applications (e.g. hyperbolic deep learning \cite{nickel2017poincare, liu2019hyperbolic}). Currently, there is \emph{no such deep learning library} for efficient high precision computations. Popular deep learning frameworks such as PyTorch \cite{paszke2019pytorch}, TensorFlow \cite{tensorflow2015-whitepaper} and Caffe \cite{jia2014caffe} only supports training with standard floats (e.g. float16, float32 and float64). Though high precision computations are enabled in Julia using multiple-digit ``BigFloat", deep learning libraries built on top of Julia such as Flux \cite{Flux.jl-2018,innes:2018}, MXNet.jl \cite{chen2015mxnet}, and TensorFlow.jl \cite{malmaud2018tensorflow} do not support training with BigFloat. Furthermore, BigFloat in Julia can \emph{only exist on CPUs}, but not on GPUs, which greatly limits its usages. 

In this work, we develop a multiple-component floating-point library, \emph{MCTensor}, for general-purpose, fast, and high-precision deep learning. We build MCTensor on top of PyTorch and it can be used in the same way as the PyTorch Tensor object. We hope MCTensor library can benefit applications in the following areas: (1) high precision training with low precision numbers. Training with low precision numbers is an emerging deep learning area on tasks like mobile computer vision \cite{tulloch2017high} and leveraging hardware accelerator like Google TPU \cite{jouppi2017datacenter}, where the deep learning model computes with low precision numbers such as float8 and float16. A large quantization error using low precision arithmetic would affect the convergence \cite{wu2018training} and may degrade the performance of the model; (2) numerical accurate and stable hyperbolic deep learning, where the non-Euclidean hyperbolic space is used in place of Euclidean space for many deep learning tasks and models due to its non-Euclidean properties. For example, graph embedding \cite{nickel2017poincare,nickel2018learning} and many developed hyperbolic networks, including hyperbolic neural networks \cite{ganea2018hyperbolic} and hyperbolic GCN \cite{chami2019hyperbolic,yu2022hyla}. However, the numerical error of computing with standard floating-point numbers in hyperbolic space is unbounded, even with float64, characterized as the ``NaN" problem \cite{sala2018representation,yu2019numerically}. A high precision computation in the hyperbolic space suggested using MCF \cite{yu2021representing} would be helpful. Our main contributions are as follows:
\begin{itemize}[nosep, leftmargin=*]
    \item We implement MCTensor in the same way as PyTorch Tensor with corresponding basic and matrix-level operations using MCF.
    \item We enable learning with MCTensor by developing the MCModule layers, MCOptimizers and etc with the same programming interface as PyTorch's counterparts.
    \item We demonstrate the performance of MCTensor for both high precision training with low precision numbers and on some hyperbolic tasks.  
\end{itemize}




 







    





\section{Related Work}
\wentao{1. high precision floats, what is bigfloat, usage of bigfloats and limitations. }

\TC{2. introduce mcf, some algorithms and guarantees proposed. requirements, orders decreasing, sum of multiple ordinary floats, advantages over bigfloat.
3.(double-double, quad-double)
existing work using mcf and how, quad-double, double-double library and other libraries, how many components, applications.}

\canal{4. applications that benefits from MCF, use MCF to achieve high precision training with low precision floats, hyperbolic space, the NaN problem, running with float64 to avoid nans, unbounded errors, why this is needed? }






High precision arithmetic were explored and implemented in many programming languages.  Bailey and Borwein \cite{bailey2015high} introduced the ongoing challenges for carrying the high precision computations, specifically in the area of mathematical physics.  Also, the software and hardware requirements for dealing with numerical difficulties were presented in light of emerging developments in computer architecture.  In this work\cite{fousse2007mpfr}, a multiple-precision binary floating-point library was presented.  The implementations were built on top of GNU MP library and the was programmed in ISO C language.    Through providing proper rounding and exceptions, the MPFR library can extend to arbitrary-precision from the IEEE 754 standard.  


A Multi-Component Floating Point Format is an unevaluated sum of multiple IEEE half/standard/double floating point numbers and can serve as a representation for higher precision floating point values.  Hida et al\cite{hida2001algorithms} showed that with a quadruple of double precision floating point values, precision of at least 212 bit of significand can be achieved.  In their work, there are also algorithms designed for carrying arithmetic operations (four basic operations and transcendental operations) in the Multi-Component format.  They also further demonstrated the extensive correctness of these algorithms.

QPyTorch\cite{zhang2019qpytorch} is one example of utilizing low precision floating point values, for instance, Half Tensor(float16), for arithmetic simulation framework based on PyTorch. Specifically, QPyTorch was designed to investigate the area of low-precision training and to aid such study through simulation without hardware overhead.  Furthermore, to reduce the impact of simulator overhead, a fused-kernel approach was deployed and thus enabling study of large-scale problems with QPyTorch.  


















\section{Methodologies}
Here we introduce some basics of our MCTensor library, built on top of PyTorch \cite{paszke2019pytorch} \footnote{Our PyTorch version is 1.11.0} that employs multi-component floating-point as its underlying tensor representation. Each MCTensor $x$ is represented as an expansion, an unevaluated sum of multiple tensors as follows:
\vspace{-2em}
\begin{equation} \label{eq:MCF}
    x = (x_0,x_1,\cdots,x_{nc-1}) = x_0 + x_1 + ... + x_{nc-1}  \vspace{-0.5em}
\end{equation}
where each $x_i$, as a component of $x$, can be a PyTorch floating-point Tensor in any precision, and $nc$ is the number of components for MCTensor $x$. It's required that all components to be ordered in a decreasing magnitude (with $x_0$ being the largest and $x_{nc-1}$ being the smallest). In this way, MCTensor allows roundoff error to be propagated to the later components and thus offers better precision compared to a standard PyTorch Tensor\footnote{unless otherwise specified, the PyTorch tensor, or ``Tensor", is referred to a PyTorch tensor with an arbitrary floating point data type in this paper}. 

We first implement basic operators \pythoninline{add}, \pythoninline{subtract}, \pythoninline{multiply}, \pythoninline{divide}, \ldots for MCTensor with MCF arithmetic and further vectorize them to matrix-level operators \pythoninline{dot, mm,} \ldots with same semantics as their PyTorch counterparts. These operators then allow us to implement higher-level \textbf{MCModule, MCOptim} as the counterpart for \pythoninline{torch.nn.Module} and \pythoninline{torch.optim} so that we can use them for any deep learning applications.

\subsection{MCTensor Object with Basic Operators}

\textbf{[MCTensor object]} A MCTensor $x$ can be abstracted as an object with specification \pythoninline{x\{fc, tensor, nc\}}. Specifically, \pythoninline{x.tensor} has $nc$ components of PyTorch tensors $x_0, x_1, \cdots , x_{nc-1}$ in the last dimension, and it has shape as $(*x_0.\text{shape}, nc)$. The $x.fc$ data term in MCTensor is a view of $x_0$, keeps track of the gradient for $x$, and if needed, serves as an approximate tensor representation of $x$.

\textbf{[Gradient]} Because a MCTensor is an unevaluated sum of Tensors, then the gradient of a function $f$ w.r.t. $x$ is $\frac{\partial f}{\partial x}  = \frac{\partial f}{\partial (x_0 + x_1 + \cdots + x_{nc-1})}= \frac{\partial f}{\partial x_i}$,
which is same as the the gradient of $f$ w.r.t. any component $x_i$.
So we only keep track of the gradient information of \pythoninline{x.fc} for a MCTensor \pythoninline{x}, which can then be computed naturally by PyTorch's auto-differentiation engine to get \pythoninline{x.fc.grad} as \pythoninline{x.grad}.

In order to develop more advanced operations of MCTensor, we develop two most basic functions first, \pythoninline{Two-Sum} and \pythoninline{Two-Prod}, whose inputs are two PyTorch Tensors and returns the result as a MCTensor with 2 components in the form of (result, error). 


\textbf{[Basic Operators]} We develop 
binary MCTensor operators to support input types (1) a MCTensor and a Tensor, (2) a Tensor and a MCTensor, and (3) a MCTensor and a MCTensor. For unary operators, we only accept a MCTensor as an input. The output for all these operators is a MCTensor with the same $nc$ as the input(s). We implement MCF algorithms for basic arithmetic, with the full list of their algorithm statements available in the appendix:
\begin{itemize}[nosep, leftmargin=*]
    \item \textbf{MCTensor:}  \pythoninline{Exp-MCN} (exp), \pythoninline{Square-MCN} (square)
    \item \textbf{MCTensor and Tensor:} \pythoninline{Grow-ExpN} (add), \pythoninline{ScalingN} (multiply), \pythoninline{DivN} (divide) 
    \item \textbf{MCTensor and MCTensor:} \pythoninline{Add-MCN} (add), \pythoninline{Div-MCN} (divide), \pythoninline{Mul-MCN} (multiply)
\end{itemize}

For example, the addition between a MCTensor and a Tensor, or \pythoninline{Grow-ExpN} (Grow Expansion with Normalization), is given in Alg.~\ref{alg:growexpn_main}.

\begin{algorithm}[H]
   \caption{\textbf{Grow-ExpN}}
   \label{alg:growexpn_main} 
\begin{algorithmic}
   \STATE {\bfseries Input:}  $nc$-MCTensor $x$, PyTorch Tensor $v$
   \STATE \textbf{initialize} $Q\leftarrow v$
   \FOR{$i=1$ to $nc$}
    \STATE $k\leftarrow nc+1-i$
    \STATE $(Q, h_k) \leftarrow \textbf{Two-Sum}(x_{k-1}, Q)$
   \ENDFOR
   \STATE $h  \leftarrow (Q,h_1,\cdots,h_{nc})$
   \STATE {\bfseries Return:} $\textbf{Simple-Renorm}(h, nc)$  
\end{algorithmic}

\end{algorithm}

\pythoninline{Grow-ExpN} takes a $nc$-component MCTensor $x$ and a normal Tensor $v$ as input. The approximated result $Q$ of this addition, is initialized to $v$. The Grow Expansion happens first in the for loop, starting with the last component of $x$, $x_{nc-1}$, to the first component $x_0$. The algorithm \pythoninline{Two-Sum} sums $Q$ and a component $x_{k-1}$, resulting in the updated approximation $Q$ and the error term $h_k$. The resulting $h$ is therefore grown into $nc+1$ components, which are naturally ordered in a decreasing manner, except that there could be some intermediate zero components. Hence, in order to meet the requirements of a MCTensor, we use \pythoninline{Simple-Renorm} algorithm to move zeros backwards and output a $nc$ components MCTensor.

We also implement two different algorithms for \pythoninline{Mul-MCN}: a fast version that uses two \pythoninline{Div-MCN} operations to perform multiplication, and a slower version that have better error bounds. However, in practice we see little difference between them, so we run all our experiments with the fast version. Details of these basic MCF algorithms are provided in the appendix.


To demonstrate how a MCTensor's precision increases with the number of components $nc$, we plot the relative numerical errors for MCTensor with different $nc$ and PyTorch Tensor. We set the data type for both MCTensors and Tensors to be Float32, and compute the errors w.r.t. the results derived by high precision Julia BigFloat with number of bits of the significand set to 3000. As can be seen from Figure~\ref{fig:mult-err}, even with $nc=2$, the relative error is orders of magnitude smaller than PyTorch Tensors. 

\begin{figure}[t]
\centering
\includegraphics[width=0.45\textwidth]{figs/julia/same/mult.png} 
\vspace{-1em}
\caption{Relative Errors for Multiplication \pythoninline{Mul-MCN(x,y)}, Compared with High Precision Julia BigFloat (3000 bits precision). Order of magnitudes for \pythoninline{x} and \pythoninline{y} are the same.} 
\label{fig:mult-err}
\vspace{-2em}
\end{figure}

\subsection{MCTensor Matrix Operators}

After defining the basic MCF arithmetic, we are able to implement commonly used matrix level operators for MCTensor including \pythoninline{AddMM-MCN} (torch.addmm), \pythoninline{Dot-MCN} (torch.dot), \pythoninline{MV-MCN} (torch.mv), \pythoninline{MM-MCN} (torch.mm), \pythoninline{BMM-MCN} (torch.bmm) and \pythoninline{Matmul-MCN} (torch.matmul). Details of them are provided in the appendix.


For matrix operators, we leverage the broadcastability and vectorization from native PyTorch operations embedded in our basic MCF operators across each $nc$, and then apply sequential error propagation. For example, for \pythoninline{MV-MCN} with input as $x$ and $v$, we first use \pythoninline{ScalingN} to compute the broadcasted product of $x$ and $v$ for each $nc$, and then we sequentially sum up the results and propagate errors with the \pythoninline{Add-MCN} operator. In this way, we can make our multiplication part (\pythoninline{ScalingN}) independent of the input size and only employ for-loop for addition part (\pythoninline{Add-MCN}) since error propagation in the addition is sequential by the algorithm. Theoretically, accurate addition of $N$ MCTensors is of order at least $O(nc\cdot\log N)$. 

MCTensor operators will be slower because of the need to propagate errors in computation, and have more memory burden than PyTorch operators because of the nature of MCF representation. In Table~\ref{tab:matrixop}, we can see a tradeoff between program speed
and precision. However, we would like to note that there is still much space to optimize these algorithms for better timing in practice. This work aims to provide ML community with the possibility to do high-precision computations for learning with MCTensor over GPUs. More details can be found in \ref{mm}.

\begin{figure*}[htbp]
\centering
{\pythonexternal{code/MCLinear.py}} \vspace{-1em}
\caption{An Example for Implementing MCLinear }
\label{fig:code_mclinear}
\end{figure*}


\begin{table*}[htbp]
 \centering
    \scalebox{0.9}{\small\begin{tabular}{lccccl} 
    \toprule
     Operators & Inputs sizes & FloatTensor & 1-MCTensor & 2-MCTensor & 3-MCTensor  \\
    \midrule
    \textbf{Dot-MCN} & $5000, 5000$ & $ 1.61 \mu\text{s}\pm3.29 \text{ns}$ 
    & $ 442 \mu\text{s}\pm5.61 \mu\text{s}$
    & $ 656 \mu\text{s}\pm1.16 \mu\text{s}$
    & $ 858 \mu\text{s}\pm12.2 \mu\text{s}$\\                
\textbf{MV-MCN} & $(5000\times 500),500$    
& $ 157 \mu\text{s}\pm4.32 \mu\text{s}$
& $ 320 \text{ms}\pm5.78 \text{ms}$
& $ 460 \text{ms}\pm10.7 \text{ms}$
& $ 580 \text{ms}\pm12.1 \text{ms}$ \\
    \textbf{Matmul-MCN} & ( 500$\times$ 200), ( 200$\times$ 50)    
    & $ 97.3 \mu\text{s}\pm1.1 \mu\text{s}$
    & $ 495 \text{ms}\pm10.8 \text{ms}$
    & $ 735 \text{ms}\pm21.7 \text{ms}$
    & $ 934 \text{ms}\pm28 \text{ms}$ \\
    \bottomrule
    \end{tabular}}
    \caption{\centering  MCTensor Matrix Operators Running Time (mean $\pm$ sd) }
    \label{tab:matrixop}
    \vspace{-1em}
\end{table*}
 


\subsection{MCModule, MCActivation, and MCOptim}

We enable learning with MCTensor by developing neural network basic modules (MCModule), activation functions (MCActivation) and optimizers (MCOptim), in the exact programming interfaces as their PyTorch counterparts in (\pythoninline{torch.nn.Module}, \pythoninline{torch.nn.functional} and \pythoninline{torch.optim}).
The semantics are identical to that of PyTorch's module and optimizer, except that we are using MCTensor arithmetic. Specifically, we give some examples:
\begin{itemize}[nosep, leftmargin=*]
    \item \textbf{MCModule:} {\pythoninline{MCLinear, MCEmbedding, MCSequential}}
    \item \textbf{MCOptim:} {\pythoninline{MC-SGD, MCAdam}}
    \item \textbf{MCActivation:} {\pythoninline{MCSoftmax, MCReLU, MC-GELU}} 
\end{itemize}
Here Fig.\ref{fig:code_mclinear} is the MCTensor implementation of linear layer, where weights and biases are represented by MCTensors. Similar to its PyTorch peer, a \pythoninline{MCLinear} takes input of size of input shape, in\_features, size of output shape, out\_features, learn with or without biases (boolean: bias) and number of component $nc$ for setting up the underlying MCTensors.  To handle the operations within MCLinear layer, we override PyTorch's \pythoninline{nn.functional.linear} to perform matrix level multiplication (Algo. \ref{alg:matmulmcn}) between a MCTensor weight and the Tensor input, then the addition of the product with bias (Algo.\ref{alg:addmcn}) if needed. The output of the MCLinear layer is a MCTensor with the same $nc$. 


\par
Since we implemented MCModule, MCActivation, and MCOptimizer to follow the same specification as their PyTorch counterparts, building a MCTensor model is the same way as one would build a PyTorch module. The only difference to the users is the need to specify the number of components $nc$. An demonstration of our MCTensor model optimization programming paradigm can be found in Fig. \ref{fig:code_mc_opt}. \vspace{-1em}
\begin{figure} [h]
\centering
\pythonexternal{code/MCOptimization.py} \vspace{-1em}
\caption{\centering MCTensor model optimization programming paradigm}
\label{fig:code_mc_opt}
\vspace{-2em}
\end{figure}
 


















 





\subsection{Error Analysis} \label{subsec:error}
Principally, one can achieve arbitrary precision using MCTensor by simply increasing the number of components, however, note that each component of a MCTensor is a standard float, which has a natural range. Take Float16 for example, the minimum representable strictly positive value is $2^{-24}\approx 5.96 \times 10^{-8}$, i.e., the smallest error that can be captured by MCFloat16 is $2^{-24}$. Hence in practice, 2-MCFloat16 is usually sufficient and similar performances are observed when more components of MCFloat16 were adopted. While simply adding an appropriate scale factor $2^{-k}$ (depending on the precision requirements) to smaller components can help capture even smaller errors, it is out of scope of this paper. 

To validate improved precision of MCTensor models, we consider the linear regression task since it is possible to obtain loss arbitrarily close to zero.
We use a single \pythoninline{MCLinear} without bias term and the Mean Squared Error (MSE) loss on fully observed synthetic data: $\mathcal{L}(W) = \text{MSELoss}(y, XW^T)$,
where $X$ is a $(10000 \times 2)$ matrix with each entries sampled from $x  \sim \mathcal{N}(-0.5, 0.5^2)$, and $y$ is calculated from $y=X{W^*}^T$ where $W^*$ is the $(10000\times1)$ target weight sampled from $w^* \sim \mathcal{N}(-0.5, 0.5^2)$. We use gradient descent with $lr=0.05$ for optimization. 


In Figure~\ref{linlossfig}, we plot the training loss curves for the model with the same structure and initialization, but with MCTensor or Tensor as data structure. The comprehensive results can be seen in Table \ref{tab:linreg}. The final train loss for 2-MCFloat16 is orders of magnitudes smaller than Float16.


\begin{figure} 
\vspace{-1em}
\centering
 \includegraphics[width=0.45\textwidth]{figs/experiments/linear.png}
\vspace{-1em}
\caption{\centering Loss Curves on Linear Regression Task.}
\label{linlossfig}
\vspace{-1em}
\end{figure}


Since by just using 2-MCFloat16, we can achieve much better precision than Float16 Tensor (HalfTensor), we run all following experiments in Float16 with $nc=$2 or 3 expect for Hyperbolic MCEmbedding. 











\section{Conclusion}
We introduce MCTensor based on PyTorch to achieve high-precision arithmetic while leveraging the benefits of heavily-optimized floating-point arithmetic. We verify its capability of high precision computations using low precision numbers and relieving the NaN problem in hyperbolic space. We hope this library could interest researchers to use MCTensor for ML applications requiring high-precision computations.

A promising future work is to design and optimize MCF algorithms to granularize the tradeoff between efficiency and precision to make MCTensor competent even for less-noisy and general tasks.
We hope this library can address the need for fast high-precision library absent for DL community and prompt DL practitioners to rethink the concept of high-precision.
\section{Appendix}
\subsection{Julia BigFloat error}
We further conduct experiments to evaluate the numerical errors of basic MCTensor arithmetic. Specifically, we first compute \pythoninline{Add-MCN(x,y)} of two random MCTensor \pythoninline{x,y} of roughly the same magnitude $m$.  \pythoninline{x,y} are sampled by first sampling two random Julia BigFloat numbers with high precision (e.g. 3000 precision) using the equation $(10-\mathcal{N}(0,1))^{m}$, then converted to their corresponding MCTensors in Float32. In order to get the \emph{exact} numerical error, we transform the MCTensor result to a Julia BigFloat number, then compute the relative error of it to the high precision addition of \pythoninline{x,y} (and not the addition of the two BigFloat numbers initially sampled) in Julia. In the same way, we compute the numerical errors of \pythoninline{Mult-MCN(x,y)} and \pythoninline{ScalingN(x,y)} with \pythoninline{y} being MCTensor in the former case and standard (PyTorch) tensor in the later case. \pythoninline{x,y} are sampled in the same way throughout the three cases. 

\begin{figure}[H]
\centering
\includegraphics[width=0.33\textwidth]{figs/julia/same/add.png}~ 
\includegraphics[width=0.33\textwidth]{figs/julia/same/scale.png}~
\includegraphics[width=0.33\textwidth]{figs/julia/same/mult.png}~
\vspace{-1em}
\caption{Relative Error of MCTensor arithmetic with different number of components, compared with high precision Julia BigFloat results (i.e. 3000 precision). Left: \pythoninline{Add-MCN(x,y)}, middle: \pythoninline{ScalingN(x,y)} and right: \pythoninline{Mult-MCN(x,y)}. Order of magnitudes for \pythoninline{x} and \pythoninline{y} are kept the same.}
\end{figure} 

For a thorough comparison, we also derive below the same numerical errors when the order of magnitudes for \pythoninline{x} varies and order of magnitudes for \pythoninline{y} is kept at 2.

\begin{figure}[H]
\centering
\includegraphics[width=0.33\textwidth]{figs/julia/varying/add.png}~ 
\includegraphics[width=0.33\textwidth]{figs/julia/varying/scale.png}~
\includegraphics[width=0.33\textwidth]{figs/julia/varying/mult.png}~
\vspace{-1em}
\caption{
Relative Error of MCTensor arithmetic with different number of components, compared with high precision Julia BigFloat results (i.e. 3000 precision). Left: \pythoninline{Add-MCN(x,y)}, middle: \pythoninline{ScalingN(x,y)} and right: \pythoninline{Mult-MCN(x,y)}. Order of magnitudes for \pythoninline{x} varies and order of magnitudes for \pythoninline{y} is kept at 2.}
\end{figure} 




\newpage
\twocolumn
\subsection{MCTensor Operators}
\subsubsection{Basic Operators}
The input of \pythoninline{Two-Sum} is two PyTorch Tensors with same precision, $a$ and $b$. Algorithm~\ref{alg:twosum} returns the sum $s = \text{fl}(a + b)$ and the error, $err(a + b)$.
\vspace{-3mm}
\input{algo/basic/TwoSum}
\vspace{-3mm}

The \pythoninline{Split} Algorithm~\ref{alg:split} takes a standard PyTorch floating point value with $p$-bit significand and splits it into its high and low parts, both with $\frac{p}{2}$-bit of significand.
\vspace{-3mm}
\input{algo/basic/Split}
\vspace{-3mm}

Based on \pythoninline{Split}, the following Algorithm~\ref{alg:twoprod} computes and returns $p = \textbf{fl}(a \times b)$ and $e = \text{err}(a \times b)$.
\vspace{-3mm}
\input{algo/basic/TwoProd}
\vspace{-3mm}

Fuse multiply-add, or FMA, is a floating-point operation that performs multiplication and addition in one step.  With proper hardware, this Algorithm~\ref{alg:twoprodfma} can speed up \pythoninline{TwoProd}.
\vspace{-3mm}
\input{algo/basic/TwoProdFMA}
\vspace{-3mm}

The constraint of decreasing magnitude and non-overlapping across $nc$ might be temporarily violated in computation, so MCTensor must be renormalized during computation.  Users can specify the target $nc$ after renormalization, $r_{nc}$, but by default we keep them the same.  The \pythoninline{Renormalize} function is this paper is a variant of the Priest's algorithm~\ref{alg:renormalize}. \cite{priest1991algorithms}.
\vspace{-3mm}
\input{algo/basic/Renorm}
\vspace{-3mm}

There is a simple and fast implementation of \pythoninline{Renormalize} as \pythoninline{Simple-Renorm}, which extracts all non-zero values from a non-renormalized MCTensor and puts together a new MCTensor.  Note that this \pythoninline{Simple-Renorm} does not have the same guarantee as \pythoninline{Renormalize}. Therefore in our \pythoninline{Mult-MCN}, we still use the original version of \pythoninline{Renormalize}.  But for most other operations, we still utilize \pythoninline{Simple-Renorm} for fast computation.
\vspace{-3mm}
\input{algo/basic/SimpRenorm}
\vspace{-3mm}

\newpage
Algorithm~\ref{alg:scalingN} describes the multiplication of a MCTensor with a PyTorch Tensor.  In our implementation, the user can specify whether the algorithm can return an expanded results with $nc+1$, or a MCTensor with same $nc$.
\vspace{-3mm}
\input{algo/basic/ScalingN}
\vspace{-3mm}

\pythoninline{Add-MCN, Div-MCN, Mult-MCN} are operators for addition, division, and multiplication of two $nc$-MCTensors. Here we have two versions of multiplication, \pythoninline{Mult-MCN} and \pythoninline{Mult-MCN-Slow}. Algorithm~\ref{alg:multmcn} is implemented by taking the inverse of the second MCTensor first, and then the first MCTensor is divided by the second MCTensor's inverse. This division would give the result of multiplication. Algorithm~\ref{alg:multmcn_slow} follows the same pattern of our definition of \pythoninline{Div-MCN} and provides better error bounds, but it is rarely used as the computational cost is too high.
\vspace{-3mm}
\input{algo/basic/AddMCN}
\vspace{-3mm}
\input{algo/basic/DivMCN}
\vspace{-3mm}
\input{algo/basic/MultMCN}
\vspace{-3mm}
\input{algo/basic/MultMCN_slow}
\vspace{-3mm}

Algorithm~\ref{alg:divN}, \pythoninline{DivN} takes input of a PyTorch Tensor and a MCTensor and computes the \pythoninline{Div-MCN} results by appending zero-value components to a PyTorch Tensor and making it MCTensor.
\vspace{-3mm}
\input{algo/basic/divn}

\newpage
The following Algorithm~\ref{alg:expn} describes the exponential function for a MCTensor.
\vspace{-3mm}
\input{algo/basic/ExpMCN}

The following Algorithm~\ref{alg:SquareN} describes the square of a MCTensor.
\vspace{-3mm}
\input{algo/basic/SquareMCN}
\vspace{-3mm}


\subsubsection{Matrix Operators}
Based on the dimensions of input MCTensor and PyTorch Tensor, \pythoninline{Dot-MCN, MV-MCN, MM-MCN, BMM-MCN and 4DMM-MCN} are implemented for calculating the matrix-level multiplication results.  All operations are identical with the PyTorch implementations.  
\vspace{-3mm}
\input{algo/matrix/DotMCN}
\vspace{-3mm}
\input{algo/matrix/MVMCN}
\vspace{-3mm}
\input{algo/matrix/MM-MCN}
\vspace{-3mm}
\input{algo/matrix/BMM-MCN}
\vspace{-3mm}
\input{algo/matrix/4DMM-MCN}
\vspace{-3mm}
The following Algorithm~\ref{alg:addmmmcn} computes the matrix multiplication of a $nc$-MCTensor with a PyTorch Tensor first, then times a constant $\alpha$.  Then the product of another $nc$-MCTensor times a constant $\beta$ is added to the former multiplication result. 
\input{algo/matrix/AddMM-MCN}
\vspace{-3mm}
\newpage
This Algorithm~\ref{alg:matmulmcn}, \pythoninline{Matmul-MCN} is the central function for handling all matrix level multiplications of one $nc$-MCTensor and one PyTorch Tensor.
\input{algo/matrix/Matmul-MCN}
\vspace{-3mm}


\onecolumn







\subsection{MCF Basic and Matrix Operators } \label{mm}
We run on a AMD Ryzen 7 5800X CPU with 64 GB memory. 

For all the basic operators, we repeat PyTorch addition for 7e3 runs and 1-, 2-, 3-MCTensor for 7 runs.

For vector-vector (\pythoninline{torch.dot}/\textbf{Dot-MCN}) product, we repeat PyTorch addition for 7e6 runs and 1-, 2-, 3-MCTensor for 7e3 runs. For matrix-vector (\pythoninline{torch.mv}/\textbf{MV-MCN}) product and batched matrix-matrix (\pythoninline{torch.bmm}/\textbf{BMM-MCN}) product, we repeat PyTorch addition for 7e3 runs and 1-, 2-, 3-MCTensor for 7 runs. For matrix-matrix (\pythoninline{torch.matmul}/\textbf{Matmul-MCN}) product and bias-matrix-matrix (\pythoninline{torch.addmm}/\textbf{AddMM-MCN}) addition-product, we repeat PyTorch addition for 7e4 runs and 1-, 2-, 3-MCTensor for 7 runs. 
 
The timing for both matrix and basic operators are given in Table \ref{tab:basicop} and Table \ref{tab:matrixop_full}.
\input{tables/basicOP}
\input{tables/matrixOP_full}

\subsection{Linear Regression Results} \label{LR}
Table \ref{tab:linreg} describes the final training loss of the linear regression task in section \ref{subsec:error}.
\input{tables/linear_train}
\newpage




 

\subsection{Logistic Regression}\label{logreg}
We apply logistic regression on a synthetic dataset and a breast cancer dataset. The synthetic dataset consists of 1,000 data points, where each data point contains two features, while the breast cancer dataset consists of 569 data points and each data point contains 30 numeric features. A single \pythoninline{MCLinear} layer with float16 and $nc$ between 1, 2 and 3 are used for logistic regression with Binary Cross Entropy loss.

$$ \mathcal{L}(W) = \pythoninline{BCELoss}(y, \text{sigmoid}(XW^T)) $$

For both datasets, we randomly split 80\% of the data for training and 20\% of the data for testing. As both datasets are small in scale, we run them in full batch with MC-SGD and SGD. For the synthetic dataset, we set learning rate to be 3e-3 and run for 4000 epochs. For the breast cancer dataset, we set learning rate to be 1e-4 and momentum as 0.9, and run for 3000 epochs. Below is the results for both dataset and figures for MCLinear results on breast cancer dataset with MC-SGD optimizer. As can be seen, even with as few as 2 components, MCTensor in float16 can match with the training loss of float32 Tensor. 
\begin{figure*}[ht] \vspace{-1em}
\centering
\begin{subfigure}
\centering 
\includegraphics[width=0.4\textwidth]{figs/experiments/logreg/logistic_cancer_test.png}  
\end{subfigure} \vspace{-1em}
\caption{\centering Testing Accuracy for MCTensor and PyTorch Logistic Regression on the Breast Cancer Dataset} 
\label{fig:logistic_cancer}
\end{figure*}
 
\begin{table*}[h]
    \centering
    \include{tables/logistic_toy}
 \vspace{-1em}
      \captionof{table}{ \centering Final Training and Testing Results for Logistic Regression on Synthetic Dataset}
\end{table*}
 
\begin{table*}[h]
    \centering
    \include{tables/logistic_cancer}
    \vspace{-1em}
      \captionof{table}{ \centering Final Training and Testing Results for Logistic Regression on Breast Cancer Dataset}
\end{table*}
 \newpage
 \subsection{MCTensor Deep Learning models} 

As a demonstration for the ease of using MCTensor to build a MCTensor deep learning model, we provide an example code for MCTensor MLP (MC-MLP) model for the multi-class classification tasks. Essentially, the only differences visible to the users are the need to set the number of components, and the need to convert MCTensor output to Tensor approximation between different layers. The need for Tensor approximation is explained below. \par

\begin{figure}[ht]
\centering
\pythonexternal{code/MLP_toy.py}
\caption{MC-MLP Code Example}
\label{fig:MLP}
\end{figure}

In MCModule, all network layers take Tensor as inputs, keep MCTensor as their weights, and through MCTensor matrix operations, produces MCTensor as output. If the MCTensor outputs need to be taken as input for the next MCModule layer, there need to be new layers that will perform multiplication on two MCTensor matrices. As can be seen from Table \ref{tab:basicop}, \textbf{ScalingN}, the multiplication between a MCTensor and a PyTorch Tensor, is more than a hundred times faster than \textbf{Mult-MCN}, the multiplication between two MCTensor. To avoid this forbidden cost of computation, we convert the unevaluated sums into an evaluated sum as shown above (\pythoninline{x.tensor.sum(-1)}). Although there might some losses of precision due to summation, the loss is marginal and we trade it for an orders of magnitude smaller execution time.

 
\subsection{MC-MLP Experiment Details}\label{MC-MLP}

Using MC-MLP, we perform multi-class classification task on the Reduced MNIST dataset and binary classification task on the Breast Cancer dataset. The training details for both tasks are shown in Table \ref{tab:mlp-b-hyper} and Table \ref{tab:mlp-m-hyper}.  
\input{tables/mlp-b-hyper}
\input{tables/mlp-m-hyper}


\newpage
Figure \ref{fig:mlp_c_sgd} shows the training loss curves for MC-MLP on breast cancer dataset with MC-SGD optimizer. 



\begin{figure}[ht]\vspace{-0.5em}
\centering
\includegraphics[width=0.4\textwidth]{figs/experiments/mlp_cancer/mlp_cancer_train.png}
\caption{\centering Training Loss Curves for MLP on Breast Cancer Dataset}
   
\label{fig:mlp_c_sgd}
\end{figure} 


Figure \ref{fig:mlp_m_sgd} shows the training loss curves for MC-MLP on reduced MNIST dataset with MC-SGD optimizer. The hyperparameters used in training are shown in Table \ref{tab:mlp-m-hyper}, and the final results for testing accuracy is shown in Table \ref{tab:mlp-m}.


\begin{figure}[!ht]
\centering
\includegraphics[width=0.5\textwidth]{figs/experiments/mlp_mnist/MNIST-MLP-train.png}
   
\vspace{-1em}
\caption{Training Loss Curves for MLP on Reduced MNIST Dataset}
\label{fig:mlp_m_sgd}
\end{figure} 

\input{tables/mlp-m}
















