algorithm
-------------

Since time does not exist *a priori* in the type IIB matrix model, solving the classical equation of motion (\[eq: EOM\]) requires some method other than just solving differential equations unlike in ordinary classical dynamics. Here we define a “potential” $$V = {\mbox{Tr}}{\left( {{\left[ A^M, {{\left[ A_M, A_0 \right]}} \right]}} - \xi A_0 \right)}^2 
+ {\mbox{Tr}}{\left( {{\left[ A^M, {{\left[ A_M, A_I \right]}} \right]}} - \zeta A_I \right)}^2 \ ,$$ which takes the minimum $V=0$ if and only if $A_M$ solves the classical equation of motion (\[eq: EOM\]). By minimizing $V$ using the gradient descent method starting from random configurations, we can generate classical solutions. Within our ansatz , the constraints reduce to $$\begin{aligned}
\label{eq:normalization}
\frac{1}{N_{\rm X}}\mathrm{Tr} {\left( X_0 \right)}^{2} = \kappa \ , 
\quad \frac{1}{N_{\rm Y}}\mathrm{Tr}M^{2} & =1 \ , \\
\frac{1}{N_{\rm X}}\sum_{i=1}^{3}
\mathrm{Tr}{\left( X_i \right)}^{2}+\frac{1}{N_{\rm Y}}
\sum_{a=4}^{9}\mathrm{Tr}{\left( Y_a \right)}^{2} & =1 \ ,\end{aligned}$$ where (\[eq:normalization\]) is chosen by using the redundancy under the rescaling $X_\mu \mapsto s X_\mu$ and $M \mapsto s^{-1} M$, which leaves $A_\mu = X_\mu \otimes M$ unaltered.

The initial configuration is prepared by generating Gaussian variables for each element of $X_\mu$, $Y_a$ and $M$ in (\[eq: ansatz of A\]). Since the equation of 