# Poseidon Hash Optimization

## Introduction

### 0. Notation & Background

#### Notation

- Let $M$ be a $3 \times 3$ MDS matrix.
- $M^{-1}$ represents the inverse of matrix $M$.
- $f: \mathbb{F}_p^3 \rightarrow \mathbb{F}_p^3$ denotes a nonlinear operation of a full round.
- $g: \mathbb{F}_p^3 \rightarrow \mathbb{F}_p^3$ signifies a nonlinear operation of a partial round.
- $\vec{s}_i \in \mathbb{F}_p^3$ refers to the state vectors of the $i$-th round.
- $\vec{r}_i \in \mathbb{F}_p^3$ indicates the round constants of the $i$-th round.
- $ARC$, $s$-$box$, and $MDS$ are functions that make up the Poseidon permutation, representing round constant addition, nonlinear operation, and MDS matrix multiplication, respectively.

#### Background

- **Sparse Matrix Factorization**: Utilized during the optimization process of partial round operations, this method factorizes a given square matrix into the product of sparse matrices (matrices largely made up of zeros). Let's define a function $D: \mathbb{F}_{p}^{3 \times 3} \rightarrow \mathbb{F}_{p}^{3 \times 3}\times \mathbb{F}_{p}^{3 \times 3}$ that takes a $3 \times 3$ matrix $M$ as input and returns two matrices $m', m''$ such that $m' \cdot m'' = M$. The function $D$ behaves as follows:

$\quad$
$$
\begin{array}{l}
\quad 0.\ \text{input } M=
\left[\begin{array}{ccc}
m_{00} & m_{01} & m_{02} \\
m_{10} & m_{11} & m_{12} \\
m_{20} & m_{21} & m_{22}
\end{array}\right]\\
\quad \\
\quad 1.\ \hat{m}=
\left[\begin{array}{cc}
m_{11} & m_{12} \\
m_{21} & m_{22}
\end{array}\right],
\hat{n}=
\left[\begin{array}{cc}
n_{11} & n_{12} \\
n_{21} & n_{22}
\end{array}\right] \text{ such that }\hat{n}=\hat{m}^{-1}\\
\quad \\
\quad 2.\ m^{\prime}=
\left[\begin{array}{ccc}
1 & 0 & 0 \\
0 & m_{11} & m_{12} \\
0 & m_{21} & m_{22}
\end{array}\right]\\
\quad \\
\quad 3.\ \vec{w}=
\left[\begin{array}{c}
m_{10} \\
m_{20}
\end{array}\right]\\
\quad \\
\quad 4.\ \vec{\hat{w}}=
\hat{n}\cdot \vec{w}\\
\quad \\
\quad 5.\ m^{\prime \prime}=
\left[\begin{array}{ccc}
m_{00} & m_{01} & m_{02} \\
\hat{w}[0] & 1 & 0 \\
\hat{w}[1] & 0 & 1
\end{array}\right]\\
\quad \\
\quad 6.\ \text{ return } m^{\prime}, m^{\prime \prime}
\end{array}\\
\quad \\
$$
This function $D$ inputs a $3\times 3$ matrix $M$ and outputs two $3\times 3$ matrices $m', m''$ such that $M = m'\cdot m''$. It's relatively straightforward to verify this.
$$
\begin{array}{l}
m'\cdot m''&=&
\left[\begin{array}{ccc}
1 & 0 & 0 \\
0 & m_{11} & m_{12} \\
0 & m_{21} & m_{22}
\end{array}\right]\cdot
\left[\begin{array}{ccc}
m_{00} & m_{01} & m_{02} \\
\hat{w}[0] & 1 & 0 \\
\hat{w}[1] & 0 & 1
\end{array}\right]\\
\quad \\
&=& \left[\begin{array}{ccc}
m_{00} & m_{01} & m_{02} \\
m_{11}\cdot \hat{w}[0]+\ m_{12} \cdot \hat{w}[1]& m_{11} & m_{12} \\
m_{21}\cdot \hat{w}[0]+\ m_{22} \cdot \hat{w}[1] & m_{21} & m_{22}
\end{array}\right]\\
\quad \\
\end{array}\\
$$
We can check $M = m'\cdot m''$ by following equations:
$$
\begin{array}{l}
\left[\begin{array}{c}
m_{11}\cdot \hat{w}[0]+\ m_{12} \cdot \hat{w}[1]\\
m_{21}\cdot \hat{w}[0]+\ m_{22} \cdot \hat{w}[1] \\
\end{array}\right]&=&
\left[\begin{array}{c}
m_{11}\ & m_{12}\\
m_{21}\ & m_{22}\\
\end{array}\right]\cdot \left[\begin{array}{c}
\hat{w}[0]\\
\hat{w}[1] \\
\end{array}\right]\\
\quad \\
&=&
\hat{m}\cdot \vec{\hat{w}}\\
\quad \\
&=&
\hat{m}\cdot \hat{n}\cdot \vec{w}\\
\quad \\
&=&
\vec{w}
\end{array}\\
\quad \\
$$

A critical point to understand is that $m'$ and $m''$ include parts of the identity matrix. This property allows us to express partial round operations using fewer constraints.
$$
\begin{array}{l}
\left[\begin{array}{ccc}
x & \cdots &\cdots\\
\end{array}\right]\cdot
\left[\begin{array}{ccc}
1 & 0 & 0 \\
0 & \ddots & \vdots \\
0 & \cdots &
\end{array}\right]=
\left[\begin{array}{c}
x \\
\vdots \\
\vdots \\
\end{array}\right]
\end{array}
$$
$\quad$

## Optimization

### 1. Full Round

A single full round consists of the operations ARC > s-box > MDS. Represented using the notation defined earlier, it looks as:
$\quad$
$$ f(\vec{s}_i+\vec{r}_i) \cdot M = \vec{s}_{i+1} $$
$\quad$
By rearranging the equation, we get:

$$(f(\vec{s}_i+\vec{r}_i) +\vec{r}_{i+1} \cdot M^{-1}) \cdot M = \vec{s}_{i+1}+\vec{r}_{i+1}$$
$\quad$
If $\vec{s}_i+\vec{r}_i$ is the output of the $i-1$ full round, the constraint for the $i$ round establishes the operational relationship between $\vec{s}_i+\vec{r}_i$ and $\vec{r}_{i+1} \cdot M^{-1}$, outputting $\vec{s}_{i+1}+\vec{r}_{i+1}$.
$\quad$
### 2. Transition from Full to Partial Round

If $\vec{s}_i+\vec{r}_i$ is the output from the initial half of the full rounds, and we continue with another full round, we can express the output using the MDS matrix $M$ as:
$\quad$
$$(f(\vec{s}_{i-1}+\vec{r}_{i-1})+\vec{r}_{i} \cdot M^{-1}) \cdot M = \vec{s}_{i+1}+\vec{r}_{i+1}$$
$\quad$
However, in the optimization method of Filecoin, a somewhat technical approach is used in the last full round to take advantage of the fact that only the first state applies a nonlinear operation in the partial round.
Assume that we have two partial rounds.
1. Operations in the Partial Round using the Sparse MDS Matrix

Let's initialize with $m'_0, m''_0 \leftarrow D(M)$ and consider $m'_1, m''_1 \leftarrow D(M \cdot m'_0)$.

The constraint for the input of the partial round in the last full round can be set as:

$$
(f(\vec{s}_{i-1}+\vec{r}_{i-1})+\vec{r}_{i}\cdot M^{-1})\cdot M\cdot m'_1 = (\vec{s}_{i}+\vec{r}_{i})\cdot {m'_1}
$$

This equation provides an understanding of the requisite operations for ARC > s-box > MDS in each partial round. Here, ARC represents the AddRoundConstant operation.

Given that $m_i'$ is of the form
$$
\left[\begin{array}{ccc}
1 & 0 & 0 \\
0 & \ddots & \vdots \\
0 & \cdots &
\end{array}\right]
$$

and since the operation $\color{red}g$ applies only to the first component, the following equation holds:

$$
(\color{red}{g(}\vec{s}_{i}+\color{green}{\vec{r}_{i}}\color{red}{)})\cdot {m'_1} = \color{red}{g(}(\vec{s}_{i}+\color{green}{\vec{r}_{i}})\cdot {m'_1}\color{red}{)}
$$

Multiplying both sides by $m''_1$ gives:

$$
\begin{array}{l}
(\color{red}{g(}\vec{s}_{i}+\color{green}{\vec{r}_{i}}\color{red}{)})\cdot {m'_1}\cdot {m''_1} = (\color{red}{g(}\vec{s}_{i}+\color{green}{\vec{r}_{i}}\color{red}{)})\cdot \color{blue}M\cdot {m'_0} = \vec{s}_{i+1}\cdot {m'_0}
\end{array}
$$

To achieve the resultant form $\vec{s}_{i+1}+\vec{r}_{i+1}$, the round constants can be modified as:

$$
\begin{array}{l}
(\color{red}{g(}\vec{s}_{i}+\color{green}{\vec{r}_{i}}\color{red}{)}+\color{green}{\vec{r}_{i+1}\cdot{M^{-1}}})\cdot {m'_1}\cdot {m''_1} = (\vec{s}_{i+1}+\vec{r}_{i+1})\cdot {m'_0}
\end{array}
$$

Repeating the same process by applying a non-linear operation to both sides and multiplying by $m''_0$, we obtain:

$$
\begin{array}{l}
(\color{red}{g(}\vec{s}_{i+1}+\color{green}{\vec{r}_{i+1}}\color{red}{)})\cdot {m'_0}\cdot {m''_0} = \vec{s}_{i+2}
\end{array}
$$

Similarly, to get the form $\vec{s}_{i+2}+\vec{r}_{i+2}$, the round constants can be modified as:

$$
\begin{array}{l}
(\color{red}{g(}\vec{s}_{i+1}+\color{green}{\vec{r}_{i+1}}\color{red}{)}+\color{green}{\vec{r}_{i+2}\cdot{M^{-1}}})\cdot \color{blue}M = \vec{s}_{i+2}+\vec{r}_{i+2}
\end{array}
$$

2. Derivation of Round Constants

The complete equation can be expressed as:

$$
\vec{s}_{i+2}+\vec{r}_{i+2} = \color{red}{g(}(\color{red}{g(}\vec{s}_{i}+\color{green}{\vec{r}_{i}}\color{red}{)}+\color{green}{\vec{r}_{i+1}\cdot{M^{-1}}}+\color{green}{\vec{r}_{i+2}\cdot M^{-2}})\cdot {m'_1}\cdot {m''_1}\color{red}{)}\cdot m''_0
$$

Given that $\color{red}g$ only operates on the first component, the round constants obtained by multiplying with the inverse MDS matrix can be precomputed. This means only the first component of the added round constant needs to be considered in the non-linear operation during the partial rounds.

Thus, defining ${\vec{r}'_{i+1}} = {\vec{r}_{i+1}\cdot{M^{-1}}}$ and $\vec{r}''_{i+2} = {\vec{r}_{i+2}\cdot M^{-2}}$, the round constants can be redefined as:

$$
\begin{array}{l}
\vec{c_i} = \vec{r_i} +
\left[\begin{array}{c}
0 \\
r'_{i+1}[1] \\
r'_{i+1}[2]
\end{array}\right] +
\left[\begin{array}{c}
0 \\
r''_{i+2}[1] \\
r''_{i+2}[2]
\end{array}\right]
\end{array}
$$

3. Partial Round Optimization: Optimized Round Constants + Sparse MDS

Thus, the constraint for the input of the partial round in the final full round is as follows:

$\begin{array}{l}
(f(\vec{s}{i-1}+\vec{r}{i-1})+\vec{c}{i}\cdot M^{-1})\cdot M\cdot m'_1 &=& (\vec{s}{i}+\vec{c}{i})\cdot {m'_1} \
\end{array}$

In the subsequent partial rounds, one doesn't need to consider the additional round constant additions for the second and third states. The focus should be on the multiplication with the sparse matrix and the addition of round constants and non-linear operations for the first state.

However, since the first component of the state involves a mixture of operations - adding round constants and non-linear computations - the round constants added to the first component when calculating $\vec{c}_i$ must be stored separately and used throughout the partial round iterations.

#### considerations:

- Utilizing the MDS matrix $M \cdot m_1'$ instead of $M$ is a consideration for executing two partial rounds. Depending on the number of rounds, there is a need for additional computations on the round constants and the matrix to be multiplied in the final full round. However, these calculations can be precomputed before execution.

- Through the optimization of round constants, the number of round constants that need to be stored, as well as the frequency of round constant operations, are reduced.

- For a generalized algorithm, please refer to Filecoin's optimization documentation.

- The techniques employing $m'$ and $m''$ can be viewed as being introduced to segregate the operations on the first component of the state vector from the second and third components.

### 3. partial round
The output from the full round to the partial round starts from $(\vec{s}_{i}+\vec{r}_{i})\cdot {m'}$. The matrix $m'$, by definition, does not change the first component of the state vector value when applying the matrix multiplication operation, and it is important to note that $g$ is a nonlinear operation applied only to the first component. After applying the nonlinear operation and multiplying the matrix $m''$, it is:
$$
g((\vec{s}_i+\vec{r}_i)\cdot m')\cdot m'' = \vec{s}_{i+1}
$$

Using the fact that $m'\cdot m'' = M$, as in the full round, we can consider the following relationship for $\vec{s}_i+\vec{r}_i, \vec{s}_{i+1}, \vec{r}_{i+1}$:
$$
(g(\vec{s}_i+\vec{r}_i)+\vec{r}_{i+1}\cdot M^{-1})\cdot M= \vec{s}_{i+1}+\vec{r}_{i+1}
$$
Through the above process, it is possible to check whether the operation for one partial round has been correctly performed, and if there are more than two partial rounds, the $M\cdot m'$ matrix used in the '2. full -> partial round' process is applied again as the input to the function $D$ as many times as the number of rounds, and only the matrix values corresponding to $m''$ are stored separately and can be extended by applying them during the partial round progression.
$\quad$
$\quad$
## Reference
[Filecoin's poseidon_spec.pdf](https://github.com/lurk-lab/neptune/blob/main/spec/poseidon_spec.pdf)
: The content above was written based on the Filecoin's document as a reference to understand the optimization technique of the Poseidon hash function. More detailed descriptions and generalized algorithms are described in the document in greater detail.