\section*{ABSTRACT}

Light emission and detection are the two fundamental attributes of optoelectronic communication systems. Until now, both functions have been demonstrated using the p-n diode which is exploited across a wide range of applications \cite{sze2021physics}. However, due to the competing dynamics of carrier injection and photocarrier collection \cite{rau2007reciprocity}, with this device light emission and detection are realized separately by switching the direction of the applied electrical bias \cite{bie2017mote,bao2020bidirectional,oh2017double}. Here we use mobile ions in halide perovskites to demonstrate light-emitting photodetection in either condition of applied electrical bias. Our device consists of a CsPbBr$_3$ microwire which is integrated with single-walled carbon nanotube thin film electrodes. The dual functionality stems from the modulation of an energetic barrier caused by the cooperative action of mobile ions with the photogenerated charge carriers at the perovskite-electrode interface. Furthermore, such complex charge dynamics also result in a novel effect: light-enhanced electroluminescence. The observed new optoelectronic phenomena in our simple lateral device design will expand the applications for mixed ionic-electronic conductors in multifunctional optoelectronic devices \cite{lee2020multifunctional}.

\hfill \break

\thispagestyle{empty}

\par The multifunctionality in optoelectronics supposes several physical phenomena involved during the device operation to serve the application demands \cite{lee2020multifunctional}. Intuitively, its possible realization is the merging of functions from separate devices combined into one architecture \cite{niu2019wireless}. However, this approach requires accurate preparation of the parts and their matching. Another approach is to optimize device operation through a controllable external stimulation. A good example is the p-n junction which can operate as a light emitter at forward bias direction and as a photodetector in the reverse bias mode \cite{bie2017mote,bao2020bidirectional,oh2017double}.
Alternatively, multifunctionality can be achieved by introducing new responsive materials to the device possessing a certain function. For instance, the exposure of p-AlGaN/n-GaN junction to electrolyte environment serving as an additional electrode enables light wavelength dependent polarity switching of photocurrent \cite{wang2021bidirectional}.

\par Along with the aforementioned approaches, there is significant interest to materials which already have structural properties affording simultaneous multiple responses. Halide perovskites are among of such materials, which in addition to their outstanding optoelectronic properties (long carrier diffusion, high carrier mobility, defect tolerance and large carrier absorption coefficient) \cite{green2014emergence, herz2017charge, stranks2013electron}, exhibit intrinsically complex mixed ionic-electronic conduction \cite{eames2015ionic}. In fact, the additive mobile ionic species in halide perovskites dynamically respond to various applied stimuli such as light \cite{kim2018large}, temperature \cite{bag2015kinetics, zhang2020defect} and electric field \cite{xiao2015giant,li2018unravelling}. Their significant impact on the device performance can be separated from electronic charge carriers \cite{kim2018large}. Indeed, the halide vacancy migration driven by applied bias modulates the local perovskite photoluminescence intensity \cite{li2018unravelling}. Furthermore, the applied electric field of 1 V/$\mu$m causes ionic redistribution in bulk perovskite leading to a switchable photovoltaic effect in metal-perovskite-conductive oxide design \cite{xiao2015giant}.

\par The sensitivity of mobile ions in halide perovskites to various stimuli is generally recognized to be the main obstacle to stable working perovskite-based common optoelectronic devices, such as solar cells (SCs), photodetectors (PDs) or light-emitting diodes (LEDs) \cite{liu2021correlations, yuan2016ion}. In these devices operating at either DC bias or built-in potential, the electric field-affected mobile ions deteriorate performance characteristics \cite{liu2021correlations}. Despite that, some emerging state-of-the-art applications avoiding light stimuli such as ionic memories and switchers take advantage of perovskite mobile ions \cite{yen2021all, zhumekenov2021stimuli}. However, the controllable ionic migration in halide perovskites upon light illumination and applying electric field has not been employed for truly multifunctional optoelectronic devices. From this point of view, unique structural properties of light-sensitive ionic semiconductors can be exploited in a rational way.

\par In this work we present a strategy to create a light-emitting photodetector operating at given electrical bias. The dual functionality (photodetection mode, light-emission mode) is achieved due to migration of perovskite ions modulating energetic barriers at the interfaces in the cooperation with the photogenerated carriers. The device is as simple as cesium-lead tribromide (CsPbBr$_3$) microwire (MW) connected to single-walled carbon nanotube (SWCNT) thin film electrodes. We demonstrate the operation of light-emitting photodetector at room temperature featuring the light-enhanced electroluminescence property .

\begin{figure}[t!]
\centering
\center{\includegraphics[width=1.0\linewidth]{Figs/Fig1.pdf}}
\caption{ | \textbf{The structure, design and operation of the CsPbBr$_3$ microwire device at 200 K}.
\\ \textbf{a} SEM image of the fabricated device: a single CsPbBr$_3$ microwire (MW) connects two SWCNT film electrodes at a distance of $\sim$ 5 $\mu m$. \textbf{b} The scheme of the device operation where the complex optical and electrical pump result in the optoelectronic dual functional response at a given polarity of applied electrical bias. \textbf{c} The schematic cross-section of MW with lead oxide (PbO) layer. The XPS depth profiling reveals thickness of PbO 4.2 $\pm$ 1.2 nm. \textbf{d,e} Current-voltage curves for the device in dark and illuminated conditions. Dark current at 200 K dramatically increases as compared to one at 160 K. The EL threshold is revealed at 4.5 V. The photocurrent is collected in the 0.6 $-$ 60 mW cm$^{-2}$ intensity range at 490 nm wavelength CW illumination. The linear fits \textbf{e} are done in range 2.25 $-$ 7 V and reflect gradual slope decrease from 4.0 $\pm$ 0.2 to 2.5 $\pm$ 0.1 with the increase in illumination. \textbf{f} Normalized EL spectra at 4.5 $-$ 7 V applied bias. \textbf{g} Schematic band diagram of the device under applied bias where the interface regions SWCNT|PbO|CsPbBr$_3$ are highlighted. The effective energetic barrier height for the carrier injection is modulated by both mobile ions $\Phi_{ion}$ and photocarriers $\Phi_{photo}$. The colors reflect the band diagram under different conditions: grey - without charge carriers, blue - with mobile ions, and red - with mobile ions and photocarriers. The corresponding red and blue horizontal lines reflect the position of the barrier for the charge carrier (electrons or holes) injection.
}
\end{figure}

\par CsPbBr$_3$ perovskite MWs are fabricated  by  temperature difference triggered growth method (for details, see Methods) \cite{wang2018temperature}. The deposition of perovskite is conducted on sapphire substrates affording the MWs directional growth (Supplementary Fig.S1) and excessive heat collection during the device operation owing to high thermal conductivity of sapphire. The typical photoluminescence spectrum collected from an ensemble of MWs is peaked at 521 nm and demonstrates full width at half maximum $\delta\lambda$ = 17.8 nm (Supplementary Fig.S1). X-ray diffraction pattern (Supplementary Fig.S2) proves high crystallinity of as-synthesized microwires possessing an orthorhombic structure (space group P$_{nma}$) \cite{saidaminov2017inorganic}. A stoichiometric atomic ratio of Cs:Pb:Br elements is checked out by the energy-dispersive x-ray analysis and equals to 1:1:3 (Supplementary Fig.S3).

\par Afterwards, SWCNT thin film synthesized by CVD method (for details, see Methods) is dry-transferred on top of the as-grown CsPbBr$_3$ MWs. Then, precise femtosecond laser ablation of the film is performed.  Generally, 1030 nm fs laser light focused at 3 $\mu$m spot with fluence below 0.1 J/cm$^2$ and 100 Hz repetition rate results in gentle removing of SWCNT film \cite{marunchenko2022single} which does not affect the surface morphology of MWs and gives an interelectrode distance of $\sim$ 5 $\mu$m (Fig.1a, Supplementary Fig.S4). The electrode function of SWCNT film is complemented by its chemical inertness which increases the operational stability of perovskite-based optoelectronic devices \cite{marunchenko2022single,aitola2017high}.

\par Eventually, the as-fabricated device is a single CsPbBr$_3$ MW short-circuiting two symmetrical laser ablated SWCNT film electrodes (Fig.1a,b, Supplementary Fig.S4). In context of semiconductor device physics, it should be referred to metal-semiconductor-metal MSM structure \cite{sze2021physics, sharma2013metal}. In general, the metal-semiconductor interface may dramatically influence the device performance. The fabrication of our device in the presence of atmosphere can lead to the formation of native metal oxide. Therefore, we perform X-ray photoelectron spectroscopy profiling (Fig.1c) which proves the formation of lead(II) oxide (PbO) of 4.2 $\pm$ 1.2 nm thickness (Supplementary Fig.S5,S6,S7). The PbO should affect carrier transport in our MSM Schottky barrier-type device by introducing an additional energetic barrier for selective charge transport \cite{sze2021physics, sharma2013metal}.

Charge carrier transport in the fabricated device is governed by coupled drift-diffusion and Poisson equations \cite{moia2019ionic,bertoluzzi2020mobile}. The latter can be written for CsPbBr$_3$ MW at coordinate $x$ as:
\begin{equation}
\frac{\partial^2{\phi}}{\partial {x^2}} = -\frac{q}{\epsilon_0\epsilon_r} (N^{+} - N^{-}),
\end{equation}
where $\phi$ is electrostatic potential, $x$ position of charges in one-dimension, $q$ is electron charge, $\epsilon_0$ and $\epsilon_r$ are the dielectric constant of vacuum and perovskite material relative permittivity, respectively. Finally, $N^{+}$ and $N^{-}$ are concentrations of positive and negative charges in perovskite. It is worth noting that the right part of \textbf{eq.1} does not reflect the nature of charge carrier type. This means that all of the present charge carriers (both ionic and electronic) govern the electric potential and electric field at both SWCNT|CsPbBr$_3$ and CsPbBr$_3$|SWCNT symmetrical interfaces.
\par For the case of metal-semiconductor (MS) interface, consider for example electron transport. The accumulation of positive carriers with concentration N$^{+}$ near MS interface region with band bending \emph{V$_{bb}$} enhances the electric field E and causes Schottky barrier $\Delta\Phi$ lowering \cite{tung2001recent,sze2021physics}:
\begin{equation}
E = \sqrt{\frac{2 q  N^+ V_{bb}}{{\epsilon_0}{\epsilon_r}}}
\end{equation}
\begin{equation}
\Delta\Phi = \sqrt[4]{\frac{2 q ^3  N^+ V_{bb}}{8 \pi^2 {\epsilon_0}^2 {\epsilon_r}^2}}
\end{equation}
This way, positive charges near MS interface N$^{+}$ lead to more efficient thermionic emission of electrons as well as enhance the probability of the tunneling through the MS interface, which leads to current increase \cite{padovani1966field,crowell1969normalized,tung2014physics}.
In our system, hybrid optoelectronic pump (Fig.1b) gives photogenerated holes and positively charged mobile ions attributing to carrier $N^{+}$. The opposite is true for the hole transport through the SWCNT|CsPbBr$_3$ interface, where similar negative charge species contribute to $N^{-}$. We further refer to the role of accumulating species ($N^{+}$,$N^{-}$) which affect carrier injection at MS interface as effective energetic barrier lowering, emphasizing mixed thermionic-field emission behaviour \cite{padovani1966field,crowell1969normalized,tung2014physics}.

\par The DC current-voltage (I-V) scans are performed in the 40 $-$ 200 K temperature range. This allows us to separate photogenerated carriers from mobile ions which possess temperature dependent transport \cite{bag2015kinetics}. The nonlinearity of I-V curves identifies the Schottky-type contacts at SWCNT|CsPbBr$_3$ interfaces (Fig.1d,e, Supplementary Fig.S8,S9). Remarkable behavior is observed when elevating temperature from 160 K to 200 K where the dark current magnitude drastically changes up to 1.8 $\cdot$ 10$^4$ times (Fig.1d, Supplementary Fig.S8). This temperature range is identified as characteristic for the mobilisation of CsPbBr$_3$ ionic species \cite{zhang2020defect}. Another evidence of the ionic motion at 200 K is the sudden change in dark current hysteresis behaviour as compared to that of observed at temperature below 160 K (Supplementary Fig.S8) \cite{lenes2011operating, yen2021all}. The mobile ions play role in effective energetic barrier lowering owing to their accumulation at SWCNT|CsPbBr$_3$ interface according to \textbf{eq.1-3}. Therefore, it is constructive to consider the effective energetic barrier modulation $\Phi_{ion}$ which gates the electronic-type current at temperatures sufficient to make the perovskite ions mobile.

\par The photocurrent invoked by continuous-wave (CW) illumination with intensity of 0.6$-$60 mW cm$^{-2}$ at wavelength 490 nm is detected as the gradual current increase (Fig.1d,1e,Supplementary Fig.S8,S9). Thus our device exhibits photodetection (PD). Notably, in the device illuminated at 200 K, not only mobile ions modulate the current via $\Phi_{ion}$, but also photogenerated carriers via $\Phi_{photo}$ do \textbf{eq.1-3}. For this case, the log($I$)/log($V$) plot in dark and light conditions is shown in (Fig. 1e) and the power law $V$ $\sim$ $I^k$ ($k$ is positive integer) is used to fit the data in the 2.25$-$7 V range. In dark conditions, the slope 4.0 $\pm$ 0.2 identifies the trap-filling regime (k > 2), whereas upon light illumination it decreases down to 2.5 $\pm$ 0.1 at 60 mW cm$^{-2}$ approaching space-charge-limited regime (k = 2) \cite{zhang2021space}.

\par In addition to photodetection, we separately reveal electroluminescence (EL) from not illuminated single CsPbBr$_3$ MW at bias above 4.5 V (Fig.1f). The appearance of EL at 200 K means the mobile ions modify the electronic component of current to a sufficient extent for making radiative recombination to be observable. The $\delta\lambda$ value of 11 nm and the peak wavelength 525 nm are achieved at 7 V and 200 K.

\begin{figure}[t!]
\centering
\center{\includegraphics[width=0.98\linewidth]{Figs/Fig2.pdf}}
\caption{ | \textbf{LEPD operation at DC at 200 K}.
\\ \textbf{a} Current on-off ratio with the integrated emission versus voltage upon 6 mW cm$^{-2}$ illumination intensity. The EL threshold separates photodetection (PD) from light-emitting photodetection (LEPD) regimes. Arrows indicate the voltage scan direction. \textbf{b,c,d} The device operation at different applied bias relevant to PD, LEPD regime and the transition between them in terms of energetic bands at interfaces SWCNT|PbO|CsPbBr$_3$. The relative concentration of charges (ions and photocarriers) at the interfaces define the relation between $\Phi_{ion}$ and $\Phi_{photo}$ gates. Electrical pump (bias) and optical pump ($h\nu$) result in LEPD causing increase in current (PD) as well as EL appearance. Blue color reflects band diagram in dark, and red upon illumination. \textbf{b} At low applied bias, the modulation $\Phi_{photo}$ is relatively high, but EL is absent. \textbf{c} Near EL threshold the operation of LEPD is defined by balanced $\Phi_{ion}$ and $\Phi_{photo}$. \textbf{d} EL of LEPD is enhanced while the relative $\Phi_{photo}$ decreases reflected in the current on-off ratio decline \textbf{a}. \textbf{e} Three regimes are separated when LEPD is illuminated: (i) EL - in dark; (ii) PL in illuminated case, when dark currrent is low; (iii) combination of PL and light-enhanced EL. \textbf{f} The novel figure of merit for LEPD shows both current and emission on-off ratio upon the illumination. The LEPD is compared to conventional p-n junction LED and PD devices.
}
\end{figure}

\par To clarify visualisation of $\Phi_{ion}$ and $\Phi_{photo}$ gating phenomena, we present the energetic bands of the device during operation (Fig.1g). The SWCNT work function (5.0 eV) as well as conduction band minimum (electron affinity) and valence band maximum for PbO (3.4 eV and 6.2 eV) and CsPbBr$_3$ (4.2 and 6.5 eV) layers are taken from previous reports \cite{shiraishi2001work,liao2016difficulty,tao2019absolute}. Here, CsPbBr$_3$ is considered to be lightly p-doped. The mobile ions together with photogenerated carriers change effective energetic barrier via $\Phi_{ion}$ and $\Phi_{photo}$. These effective enegetic barrier lowering is controlled by complex optoelectronic stimuli which result in the device multifunctionality.

\par The discovery of separate PD and EL phenomena in one device at similar bias grounds the search towards intermediate operation regime where they exist simultaneously. We introduce a term light-emitting photodetector (LEPD) for a device which enables this regime. The LEPD regime can be vividly determined by synchronizing electrical pump with a spectrometer, the information on EL at every point of the current-voltage scans is obtained (for details, see Methods). The transition between the pure PD and LEPD occurs at threshold bias for the electroluminescence (Fig.2a). Here we plot current on-off ratio and integrated light emission at 6 mW cm$^{-2}$ illumination intensity versus applied bias. The current on-off ratio is a characteristic of LEPD which separates purely ionic and mixed ionic-photo gating phenomena (Supplementary Fig.S10).

\par In LEPD regime $\Phi_{ion}$ and $\Phi_{photo}$ cooperatively affect the effective energetic barrier for carrier injection at the interfaces. However, since the gates are defined by accumulative charge $N^{+}$($N^{-}$) (eq.2-3), they compete in terms of their relative contribution to the current $I$ enhancement (Supplementary Fig.S10). We can introduce effective energetic barrier $\Phi(V)$ and consider $I$ $\sim$ $e^{-\Phi(V)}$ similar to ideal diode equation with effective energetic barrier $\Phi(V)$ \cite{sze2021physics,sze1971current,tung2014physics}. For the case of DC scans with increasing voltage (Fig. 2a), the relative concentration of photo- to ionic carriers at interfaces (Fig.2b-d) determine the relation $\Phi_{photo}$ /$\Phi_{ion}$.

\par From the point of view of band-bending, the PD, LEPD regimes and transition between them can be depicted by three diagrams showing all types of charge carriers involved to $\Phi_{ion}$ and $\Phi_{photo}$ gatings during the device operation (Fig. 2b-d) in the 1.5-7.0 V range of applied bias. At a given $\Phi_{photo}$, the $\Phi_{ion}$ increases (under bias and time) enhancing total light emission, but affecting current on-off ratio by decrease of effective energetic barrier change $\Phi(V)$ under $\Phi_{photo}$. At lower voltages PD regime is present with maximum current on-off ratio (Fig. 2a,b). Here EL is absent due to subbandgap voltages and low $\Phi_{ion}$ (Fig. 2e). In contrast, at high voltages relative $\Phi_{photo}$ gate is reduced leading to current on-off ratio decrease down to 1.8 (Fig. 2a,2d). This high-voltage LEPD regime yields EL enhanced by strong barrier modulation because of perovskite mobilie ions $\Phi_{ion}$ (Fig. 2g). Eventually, there is an intermediate regime, in between, where $\Phi_{ion}$ and $\Phi_{photo}$ gates have comparable values $\Phi_{photo}\sim\Phi_{ion}$ so that current on-off ratio is mediate, and EL threshold is overcome (Fig. 2a,2c).

\par A remarkable feature of the device operating in LEPD regime is improvement of EL intensity in the presence of incident light (Fig. 2e). The feature manifests in the reduction of turn-on voltage from 4.5 V  to 2.75 V (Fig. 2e). This means that in LEPD regime, both stimuli, namely applied voltage and incident light simultaneously enhance radiative recombination via current increase invoked by $\Phi_{ion}$ and $\Phi_{photo}$. At first sight, it is reasonable to assume that photoexcited PL makes a major contribution to the established enhancement since it is not straightforward to separate PL and EL signals (Supplementary Fig.S11). However, it should be noted that PL intensity upon 6 mW cm$^{-2}$ illumination is an order of magnitude lower than EL intensity at turn-on voltage. Moreover, a ten-fold increase in intensity of the incident light up to 60 mW cm$^{-2}$ does not affect the enhancement. Thus, we conclude that despite both mechanisms may occur, the light-enhanced EL is supposed to be main contributor to the emission improvement.

\par The emergence of new-type devices requires a reliable criteria for the comparison of their performance. As one of such criteria, we suggest figure of merit (FoM) showing emission on-off ratio versus current on-off ratio upon illumination and above EL threshold (Fig. 2f). The proposed FoM reflects the difference between the operational ranges for conventional devices (LED,PD) and LEPD. According to this, a high performance LEPD should demonstrate maximum on-off ratios at minimum applied stimuli (bias, illumination). In (Fig. 2f) one can see that on-off ratios decrease with the increase in applied voltage. Such behavior is explained by recombination kinetics in halide perovskites. In detail, when numerous trap states are present at the interfaces, they cause the accumulation of photogenerated charge carriers and contribute to $\Phi_{photo}$. At low bias, this contribution is the highest because of unfilled trap states. The latter is confirmed by the slope of log-log plot of EL intensity versus current which equals 2 and corresponds to trap-filling dominant recombination regime (Supplementary Fig.S12).

\par The temporary stable response in LEPD requires stationary $\Phi(V)_{ion}$, however dynamic accumulation of ions at the interfaces results in time-dependent $\Phi(V,t)_{ion}$ and is a major source of device instability at DC bias. A similar statement is true for other perovskite-based p-n junction-type devices \cite{eames2015ionic,liu2017lab,yuan2016ion}. Taking this into account, after the discovery of the dual functionality in the DC-biased device which additionally undergoes overheating because of Joule losses, we consider pulsed biasing affording quasi-stationary $<\Phi(V)_{ion}>$ for the device operation at room temperature. Furthermore, such AC biasing helps to resolve LEPD regime in time (Fig.3a).

\begin{figure}[t!]
\centering
\center{\includegraphics[width=0.99\linewidth]{Figs/Fig3.pdf}}
\caption{ | \textbf{LEPD operation at AC at room temperature}.
\\ \textbf{a} Scheme of the experimental setup which is used during device operation. At room temperature AC regime is provided by wave generator of varied voltage amplitude A, pulse frequency $\nu$ and Duty cycle (Dc). The CW laser source is chopped with corresponding optical frequency $\nu$. The optoelectronic signal is measured by oscilloscope and spectrometer. \textbf{b} Equivalent scheme of CsPbBr$_3$ MW presented as two Schottky diodes reflecting the interfaces modulated by time-dependent $\Phi(V,t)_{ion}$ and $\Phi(V,t)_{photo}$. The bulk of CsPbBr$_3$ perovskite provides the interfaces with mobile ions and photocarriers. \textbf{c} The transition to LEPD regime under AC rectanglur pulses is visualised at $n_{th}$ and $n+1_{th}$ pulse.
Here $\Phi(V)_{ion}$ dynamically increases in time with Dc at fixed pulse frequency. At the end of the pulse $\Phi(V)_{ion}$ relaxation is plotted with dotted curve. The difference $\Delta\Phi_{ion}(Dc)$ is Dc dependent change of $\Phi_{ion}$ between the pulses. \textbf{d-f} The 10 Hz, 6V and Dc 30 $\%$ chosen as optimal for LEPD demonstration at AC \textbf{f}. \textbf{d} $\Phi(V)_{ion}$ is changed during the pulse length leading to EL threshold at 30$\%$ Dc. \textbf{e} The dependence of rate of $\Phi(V)_{ion}$ versus V is demonstrated and the EL threshold appears at 5 V pulse amplitude. \textbf{f} The LEPD operation upon 6 mW cm$^{-2}$ illumination. The current on-off ratio 1.5 with peak emission on-off ratio of 8 (inset image) is achieved. The MW is visualised with its emission area position and the bias polarity (inset image).
}
\end{figure}

\par The LEPD MSM device can be considered as two back to back Schottky diodes which can be modulated and are connected through the perovskite resistive element (Fig.3b). Each of the Schottky diode is the interface between metal (SWCNT film) and semiconductor (CsPbBr$_3$) including PbO layer \cite{sze2021physics, tung2014physics}. Under electrical and optical stimuli, the bulk CsPbBr$_3$ provides the Schottky diodes with photo- and ionic carriers and modulates the total current via $\Phi_{photo}$ and $\Phi_{ion}$, respectively. The rate of the modulation is then defined by drift of these carriers from bulk towards the interface \cite{wang2019kinetic,bernards2007steady}.

\par For the AC biasing, periodic rectangular pulse train with zero bias between the pulses is employed (Fig. 3a,3c),(for details, see Methods). The stimulus is defined by three parameters: voltage amplitude $V$, pulse repetition frequency $\nu$ = $1/T$ and pulse width $\tau$. We schematically illustrate biasing conditions impact on the $\Phi(V)_{ion}$ (Fig. 3c). During the $n_{th}$ pulse of width $\tau$ and voltage amplitude $V$, ions modulate $\Phi(V)_{ion}$. In this case, the time response of LEPD is defined by the speed of reaching $\Phi(V)_{ion}$ sufficient for EL governed by a single-pulse width $\tau$. We identify this operation mode as a single-pulse one. In general, this value may not exceed the limit for LEPD at this pulse width $\tau$ (Fig. 3c). At the end of the pulse, where voltage is not applied, ions relax to their initial position during the timescale of $T - \tau$ which is reflected in $\Phi(V)_{ion}$ decrease (Fig. 3c). At the start of the next ${n+1}_{th}$ pulse, ions may be not completely relaxed during T - $\tau$ interval leading to small $\Delta\Phi_{ion}$ as compared to previous pulse. The ${n+1}_{th}$ pulse reaches the limit for EL threshold and, thus, LEPD regime. Here, LEPD works in accumulative mode and its time response is defined by ${n+1}_{th}$ pulse. We stress that entering LEPD regime can be realized in wide range of pulse repetition frequency, but for a given frequency, the Duty cycle (Dc) = $\nu\tau$ should be chosen according to ionic relaxation dynamics.

\par In (Fig. 3d) the dependence of the current flowing through the series shunt resistance on Dc at frequency 10 Hz and amplitude 6 V is shown. Starting from 2 {$\%$} the Dc is increased up to 30 {$\%$} until the observation of EL. The speed of the ionic modulation $\frac{d\Phi(V)_{ion}}{dt}$ at 30 $\%$ Dc increases with the applied voltage amplitude (Fig. 3e). The EL threshold is achieved at 5 V. In principle, EL is also observed at higher frequencies (e.g. 1 kHz and 100 kHz) relevant to accumulative mode of LEPD (Supplementary Fig.S13). Finally, for the set of AC biasing parameters (10 Hz, 6 V, Dc $\%$) the LEPD operating upon 6 mW cm$^{-2}$ illumination intensity is demonstrated. The photodetection is observed with on-off ratio 1.5 together with light-enhanced electroluminescence with peak emission on-off of 8 (Fig. 3f). This result states the principal ability of both current (PD) and emission (Light-enhanced EL) enhancement under the optical stimulus in one bias direction.
\\
\par  The observed phenomena of the light-emitting photodetection and light-enhanced electroluminescence in the simple MSM structures of CsPbBr$_3$ MWs underpins the unique structural complexity governed by mixed ionic-electronic conduction behaviour in the family of halide-perovskites. We expect that the cooperative action of charge carrier species (ions, photogenerated carriers) which leads to simultaneous light detection and emission will expand the applications towards multifuncitonal wearable and implantable devices needed for healthcare monitoring (e.g. pulse oximeters, smart contact lenses, optogenetics systems, neural monitoring) as well as will continue to develop in the advanced neuromorphic computing applications (which requires the memory function provided by ions) complemented by photodetection and optical switching (light-enhaced electroluminescence). The MSM halide-perovskite based devices will contribute to exploration of ionic-switching devices (e.g. ionic memories, ionic transistors), while the accurate analysis of this structure will help the researchers in the analysis of traditional optoelectronic-based devices dealing with materials with mixed ionic-electronic conduction behaviour.

\section*{Methods}

\textbf{Synthesis of CsPbBr$_3$ MWs}
\par We utilized temperature difference triggered growth method for the CsPbBr$_3$ MWs synthesis \cite{wang2018temperature}. The furnace (PZ 28-3TD High Temperature Titanium Hotplate and Program Regler PR5-3T) was used to controll the temperature during MWs growth. The CsPbBr$_3$ material was sublimated from the source substrate to the target substrate. On the glass source substrate we synthesised dense CsPbBr$_3$ microcystals according to the protocol \cite{pushkarev2018few}. The target substrate was a cleaned (5 min water and afterwards 5 min isopropanol places in ultrasonicat bath) sapphire. Two substrates were separated with air gap of 1 cm. The temperature of both substrates was controlled by the furnace temperature. The synthesis started at a furnace temperature temperature 350$^o$C. Then, the temperature rose up to 520$^o$C for 10 min and remained at this temperature for 10 min. As a result, the sublimated CsPbBr$_3$ material assembled on the sapphire substrates according to sapphire crystallographic planes.
\textbf{Device fabrication}
\par Aerosol chemical vapour deposition was implemented to synthesize SWCNT \cite{tian2011controlled}. The as-fabricated SWCNT thin films contained both metallic and semiconducting SWCNT. The film thickness of 18 $\pm$ 5 nm was estimated by combined optical absorbance and ellipsometry measurements. Later, the films were cut with a blade into stripes of width 5 mm and length 7 mm. The stripes of single-walled carbon nanotube thin films were dry-transferred on the as-synthesized CsPbBr$_3$ microwire  sapphire substrates. Afterwards, 10 $\mu$l of isopropanol (IPA) was dropped on top of SWCNT and was annealed at 80 C$^{o}$ for 5 minutes. Then the ablation was performed with light conversion pharos femtosecond laser (200 fs pulse duration, wavelength 1030 nm). The light was focused through 50x lens (NIR Mitutoyo, NA = 0.65) while the xy position was altered with a programmable piezocontroller (Standa). In such geometry, where SWCNT stripe is located above the CsPbBr$_3$ microwire, the ablation fluence should exceed the threshold fluence for SWCNT material removal while at the same time should be less then two-photon ablation of CsPbBr$_3$ material. Hence, 1) pulse repetition rate was decreased down to 100 Hz and 2) the fluence level was controlled to be below 0.1 J/cm$^2$. The laser cut was performed in two steps. Firstly, the CsPbBr$_3$ microwire was chosen and the aforementioned ablation regime for SWCNT cut above the microwire was performed. Secondly, the fluence level was increased to ablate whole SWCNT stripe and to avoid other CsPbBr$_3$ microwire interconnects in the interelectrode space.
\textbf{Optoelectronics measurements}
\par Optical measurements were performed using photoluminescence (PL) and electroluminescence (EL) spectroscopy in our custom-built experimental setup. For PL measurements the laser was focused on the samples using objective lens, the laser spot diameter was  30 $\mu m$. The excited signal (PL,EL) was collected by the same microscope objective. To excite samples we used 490 nm CW laser, to cut off excitation signal we used 500 nm long pass filter. The excited signal was measured by projecting the image plane of the objective onto the entrance slit of the imaging spectrometer (Princeton SP 2550), and detecting the light with a CCD detector (PyLoN 400BR eXcelon). The sample was mounted in an ultra-low-vibration closed-cycle helium cryostat (Advanced Research Systems, DMX-20-OM), sample holder was augmented with electrical pins for optoelectronics measurements. The optical measurements were performed on the same optical setup for AC and DC electrical measurements.

\par We used two different electrical schemes for dynamic (AC) and current-voltage curves (DC) measurements, which connected to the sample via the same port in cryostat. To measure current-voltage curves in DC regime we used Keithley 2401 Source meter. Instrument applied voltage of bias of 20 ms pulse, and the obtained current signal which was averaged through 1 second period. The scan rate was fixed to be 0.25 V per second.
Optical signal collected the same way and the illumination intensity of 490 nm CW laser was tuned through the optical filter ranged from 0.06 to 60 mW cm$^{-2}$. To electrically excite our sample in AC, we used wave generator (Keysight 33600A), which allowed us to vary duty cycle, voltage amplitude and pulse repetition frequency independently. The current signal was collected through the shunt resistance $R_{shunt}$ = 8.2 kOm on the Oscilloscope (Keysight DSOX6004A). In case of 10 Hz measurements, light illumination was chopped accordingly pulse repetition frequency, so that illumination turned-on on the scale of applied voltage pulse. The electrical measurements coordinated in time with the spectrometer, so that the device response (e.g. EL,PD) was measured simultaneously.


\section{Conclusions}
\label{sec:conclusions}

This paper presents an algorithm for generating and solving shape constraints for Deep Learning
graphs with dynamic shapes. This algorithm allows either fully static or partially static shapes to
be discovered within a program, based on shape information provided by the programmer, and the
structure of the graph. These shapes can include arithmetic expressions for their individual
dimensions, and can also vary in rank. This is done using a fully automated shape constraint solver.

\pagebreak

\section{Shape Constraint Generation}
\label{sec:constraint-generation}

Our approach extends Hindley-Milner type inference~\cite{Damas1982TypeInference}, using standard
type unification with an additional rule which generates shape constraints. This rule is described
in Section~\ref{sec:constraint-generation-rules}. From these generated shape constraints, a shape
solver (based on unification) computes concrete shapes where possible from these shape
constraints. This part of the algorithm is covered in detail in Section~\ref{sec:shape-solver}.

Given a program, which is functional in style and may be annotated with shapes by the programmer,
Hindley-Milner type inference with standard unification is run over the program in order to infer
the types of all functions, values and let bindings. During this process, the additional rule in the
unification algorithm generates a set of shape constraints for the program. A shape solver is run
over these constraints in order to solve them. Any shape variable assignments that the shape solver
discovers are then replaced in the types of the program. This produces programs with concrete shape
information where enough information is present to deduce them.

\subsection{Constraint Generation Rules}
\label{sec:constraint-generation-rules}

In order to infer information about shapes within a program, a set of shape constraints first needs
to be generated for a program. This is done as part of the standard Hindley Milner style type
inference and unification algorithms.

Unification is augmented to generate not only a type substitution for a pair of types, but also a
set of type constraints. Whenever two tensor types involving shapes are unified, shape constraints
are generated and added to the set. After unification is complete these shape constraint sets are
solved using the shape solving algorithm presented in section~\ref{sec:shape-solver}. This shape
constraint generation and solving is done incrementally as part of the type inference algorithm.

The additional unification rules for tensor types are given in
Figure~\ref{fig:unification-tensor-rules}. $G$ denotes the set of types to be unified (initially $G =
\{ \tau_1 = \tau_2 \}$ where $\tau_1$ and $\tau_2$ are the two types being unified. $S$ denotes the
set of shape constraints, initially $S = \emptyset$. If any rules match $\langle G, S \rangle$ and
generate a $\text{FAIL}$ then unification fails. Unification finishes when none of the rules match
anything in $\langle G, S \rangle$. At this point, $G$ will contain a substitution from type
variables to types, and $S$ will contain the set of shape constraints.

The first rule in Figure~\ref{fig:unification-tensor-rules} just reorders the type so that type
variables (denoted $a$) appear on the left hand side. The second rule generates a new shape
constraint given two tensor types with shapes $s_1$ and $s_2$. The third rule causes unification to
fail if a free type variable appears on both the left hand side (denoted $a$) or anywhere in the
type on the left hand side (denoted $\tau$).

\begin{figure}
  \begin{align*}
    \langle G \cup \{ \tau s = a \}, S \rangle
    & \Rightarrow \langle G \cup \{ a = \tau s \}, S \rangle \\
    \langle G \cup \{ \tau_1 s_1 = \tau_2 s_2 \}, S \rangle
    & \Rightarrow \langle G \cup \{ \tau_1 = \tau_2 \}, S \cup \{ s_1 = s_2 \} \rangle \\
    \langle G \cup \{ a = \tau s \}, S \rangle
    & \Rightarrow \text{FAIL} \enspace \text{if} \enspace a \in \text{ftv}(\tau)
  \end{align*}
  \caption{Unification rules for tensor types to generate shape constraints}
  \label{fig:unification-tensor-rules}
\end{figure}

\section{Example Graphs}
\label{sec:example-graphs}

This section presents example graphs, and how shape inference can determine the shapes of
their inputs and outputs.

\subsection{Attention}

A commonly used computation in machine learning models is attention \cite{Vaswani2017Attention}.
This can be written as the following mathematical function:

\begin{equation*}
\text{attention}(Q,K,V) = \text{softmax}(\frac{QK^T}{\sqrt{d_k}}V)
\end{equation*}

This can then be implemented in Axon as follows:

\begin{lstlisting}
softmax = fn (e : t [m,n]) {
  z = exp(e);
  s = map(
    fn (a) {
      reduce (+, 0f, a)
    },
    transpose(e)
  );
  map (fn(a) { a / s }, z)
};

attention = fn(Q, K, V) {
  matmul(softmax(matmul(Q,transpose(K))), V)
};
\end{lstlisting}

In this program, all of the types and shapes are implicit and unconstrained, except for the input to
the softmax function, which has type $t \term{[} m \term{,} n \term{]}$. This type annotation
constrains the input to be of any element types, with a shape which is two dimensional, with
dimensions $m$ and $n$ respectively.

Axon includes rules for generating shape constraints for all of the built-in functions used in this
program, such as matmul, transpose and map as described in
Section~\ref{sec:example-shape-constraints}. These rules generate a set of constraints, which when
solved give the following type and shape signature for the attention function:
\begin{equation*}
  t \term{[} a \term{,} b \term{]} \term{,} \enspace t \term{[} c \term{,} b \term{]}
  \term{,} \enspace t \term{[} c \term{,} d \term{]} \enspace \rightarrow \enspace t \term{[} a
    \term{,} d \term{]}
\end{equation*}

\subsection{Loops and Scans}

Loops and scans, as used in LSTM models \cite{LSTM}, can be written in Axon as follows:

\begin{lstlisting}
h_out, y = loop(h_in, cell, x)
\end{lstlisting}

where \lstinline|cell| is a function that takes the current state and outputs the next state,
\lstinline|x| is the input tensor, \lstinline|h_in| is the initial hidden state, \lstinline|h_out|
is the final hidden state and \lstinline|y| is the output tensor.

This loop construct can be used to write an LSTM as follows, by simply specifying the appropriate
cell function for the LSTM:

\begin{lstlisting}
cell = fn (state : (t s1, t s2), y : t s3) {
  h1, c1 = state;
  hnext, cnext = lstmStep(
    y, h1, c1,
    Wi, Wf, Wo, Wc,
    Ri, Rf, Ro, Rc,
    Bi, Bf, Bo, Bc);
  (hnext, cnext), cnext
};
loop((h, c), cell, x)
\end{lstlisting}

As with the attention example, most of this program is left with unconstrained types and shapes. The
only requirement is that the inputs to the cell function are a tuple of two tensors, and a
tensor. In this example the shape solver will determine that the output of the loop has the same
shape as the initial input tensor x, and given the shape of x, constrains the cell function to have
the correct shapes. Namely, that $s1$, $s2$ and $s3$ should be equal to the shape of x.

\subsection{Recurrent Neural Networks}

Bidirectional recurrent neural networks can be written in Axon as follows:

\begin{lstlisting}
bidir = fn (a, b, c, d, s) {
    _, x = loop(a, b, s);
    _, y = loop(c, d, reverse(s));
    concat(x, reverse(y))
};
Y = bidir(s0, A, s0', A', X);
\end{lstlisting}

This uses the same loop construct as the LSTM example previously.

Similarly, residual recurrent neural networks can be written in Axon as follows:

\begin{lstlisting}
residual = fn (a, b, s) {
  _, r = loop(a, b, s);
  r + s
};
X1 = residual(..., LSTM1, X0);
X2 = residual(..., LSTM2, X1);
\end{lstlisting}

\section{Inference Examples}
\label{sec:examples}

This section presents a set of example programs, and how shape inference can determine the shapes of
their inputs and outputs, based on both shapes provided by the programmer and the structure of the
program itself.

\subsection{Basic Symbolic Constraints}

The following Axon function performs matrix multiplication on a pair of square matrices of size
$n$. The inputs have shape $[n,n]$ and the return type is given name $r$. This example demonstrates
how the output shape is inferred from the input shape, and the structure of the program:
\begin{lstlisting}
fn (x : f32 [n,n], y : f32 [n,n]) : r {
  matmul(x, y)
}
\end{lstlisting}

The program uses the generalized matrix multiplication operator $\term{matmul}$, whose type
specification is as follows, where $\tau$ denotes the element type of the tensor:
\begin{equation*}
  (\tau \enspace s \term{@} \term{[} d_1 \term{,} d_2 \term{]},
  \enspace \tau \enspace s \term{@} \term{[} d_2 \term{,} d_3 \term{]}) \rightarrow \enspace
  \tau \enspace s \term{@} \term{[} d_1 \term{,} d_3 \term{]}
\end{equation*}

The set of shape constraints, generated by type checking this program, are:
\begin{align*}
  s \term{@} \term{[} d_1 \term{,} d_2 \term{]} & = \term{[} n \term{,} n \term{]} \\
  s \term{@} \term{[} d_2 \term{,} d_3 \term{]} & = \term{[} n \term{,} n \term{]} \\
  s \term{@} \term{[} d_1 \term{,} d_3 \term{]} & = r
\end{align*}

The first two constraints match the empty shapes rule, which assigns $s$ to be the empty shape. The
constraints are now:
\begin{equation*}
\begin{aligned}[c]
  \term{[} d_1 \term{,} d_2 \term{]} & = \term{[} n \term{,} n \term{]} \\
  \term{[} d_2 \term{,} d_3 \term{]} & = \term{[} n \term{,} n \term{]}
\end{aligned}
\quad
\begin{aligned}[c]
  s \term{@} \term{[} d_1 \term{,} d_3 \term{]} & = r \\
  s & = \term{[} \term{]}
\end{aligned}
\end{equation*}

The unpacking rule for shapes of known rank now applies to the first two constraints, giving:
\begin{equation*}
\begin{aligned}[c]
  d_1 & = n \\
  d_2 & = n \\
  d_3 & = n
\end{aligned}
\quad
\begin{aligned}[c]
  s \term{@} \term{[} d_1 \term{,} d_3 \term{]} & = r \\
  s & = \term{[} \term{]}
\end{aligned}
\end{equation*}

The reordering rule matches the fourth constraint, giving:
\begin{equation*}
\begin{aligned}[c]
  d_1 & = n \\
  d_2 & = n \\
  d_3 & = n
\end{aligned}
\quad
\begin{aligned}[c]
  r & = s \term{@} \term{[} d_1 \term{,} d_3 \term{]} \\
  s & = \term{[}\term{]}
\end{aligned}
\end{equation*}

The basic rule now matches constraint five, replacing $s$ with the empty shape in constraint four:
\begin{equation*}
\begin{aligned}[c]
  d_1 & = n \\
  d_2 & = n \\
  d_3 & = n
\end{aligned}
\quad
\begin{aligned}[c]
  r & = \term{[} \term{]} \term{@} \term{[} d_1 \term{,} d_3 \term{]} \\
  s & = \term{[} \term{]}
\end{aligned}
\end{equation*}

The basic rule also matches constraints one and three, replacing $d_1$ and $d_3$ in constraint four:
\begin{equation*}
\begin{aligned}[c]
  d_1 & = n \\
  d_2 & = n \\
  d_3 & = n
\end{aligned}
\quad
\begin{aligned}[c]
  r & = \term{[} \term{]} \term{@} \term{[} n \term{,} n \term{]} \\
  s & = \term{[} \term{]}
\end{aligned}
\end{equation*}

Shape simplification then removes the empty shape append in constraint four:
\begin{equation*}
\begin{aligned}[c]
  d_1 & = n \\
  d_2 & = n \\
  d_3 & = n
\end{aligned}
\quad
\begin{aligned}[c]
  r & = \term{[} n\term{,} n \term{]} \\
  s & = \term{[} \term{]}
\end{aligned}
\end{equation*}

Shape inference is now complete, as no rules match any of the constraints. The return shape has been
correctly inferred as $\term{[}n\term{,}n\term{]}$. Furthermore, the remaining constraints form a
shape substitution, mapping variables to shapes. The shapes for the arguments and return from the
$\term{matmul}$ call can then be found by applying this substitution to the type of
$\term{matmul}$):
\begin{equation*}
  (\tau \term{[} n \term{,} n \term{]}, \tau \term{[} n \term{,} n \term{]}) \rightarrow \tau \term{[} n \term{,} n \term{]}
\end{equation*}

\subsection{Arithmetic Expressions for Shapes}

The following program performs a 2d convolution of an 8x8 filter over an input of unknown shape
$i$. The result shape is known. This example demonstrates that shape inference can infer the shape
of the input from the shape of the output.
\begin{lstlisting}
fn (x : t i,
    f : t [4,8,8,8]
): t [4,8,1024,256] {
  conv(x, f)
}
\end{lstlisting}

The set of type constraints generated by type checking is:
\begin{equation*}
\begin{aligned}[c]
  \term{[} n \term{,} c \term{,} h \term{,} w \term{]} & = i \\
  \term{[} k \term{,} c \term{,} r \term{,} s \term{]} & = \term{[}4\term{,}8\term{,}8\term{,}8\term{]} \\
  \term{[} n \term{,} c \term{,} (\term{+} \enspace 1 \enspace (\term{-} \enspace h \enspace r))\term{,}
            (\term{+} \enspace 1 \enspace (\term{-} \enspace w \enspace s))\term{]}
  & = \term{[}4\term{,}8\term{,}1024\term{,}256\term{]}
\end{aligned}
\end{equation*}

The second and third constraints matches the unpacking rule for shapes of known rank, which gives
the following set of constraints:
\begin{equation*}
\begin{aligned}[c]
  \term{[} n \term{,} c \term{,} h \term{,}  w \term{]} & = i \\
  n & = 4 \\
  k & = 4 \\
  c & = 8
\end{aligned}
\quad
\begin{aligned}[c]
  r & = 8 \\
  s & = 8 \\
  (\term{+} \enspace 1 \enspace (\term{-} \enspace h \enspace r)) & = 1024 \\
  (\term{+} \enspace 1 \enspace (\term{-} \enspace w \enspace s)) & = 256
\end{aligned}
\end{equation*}

The basic rule for variable replacement then produces:
\begin{equation*}
\begin{aligned}[c]
  \term{[} 4 \term{,} 8 \term{,} h \term{,}  w \term{]} & = i \\
  n & = 4 \\
  k & = 4 \\
  c & = 8
\end{aligned}
\quad
\begin{aligned}[c]
  r & = 8 \\
  s & = 8 \\
  (\term{+} \enspace 1 \enspace (\term{-} \enspace h \enspace 8)) & = 1024 \\
  (\term{+} \enspace 1 \enspace (\term{-} \enspace w \enspace 8)) & = 256
\end{aligned}
\end{equation*}

Reordering reorders the first constraint:
\begin{equation*}
\begin{aligned}[c]
  i & = \term{[} 4 \term{,} 8 \term{,} h \term{,} w \term{]} \\
  n & = 4 \\
  k & = 4 \\
  c & = 8
\end{aligned}
\quad
\begin{aligned}[c]
  r & = 8 \\
  s & = 8 \\
  (\term{+} \enspace 1 \enspace (\term{-} \enspace h \enspace 8)) & = 1024 \\
  (\term{+} \enspace 1 \enspace (\term{-} \enspace w \enspace 8)) & = 256
\end{aligned}
\end{equation*}

The simplification algorithm applied to the arithmetic expressions reorders them into canonical
form:
\begin{equation*}
\begin{aligned}[c]
  i & = \term{[} 4 \term{,} 8 \term{,} h \term{,} w \term{]} \\
  n & = 4 \\
  k & = 4 \\
  c & = 8
\end{aligned}
\quad
\begin{aligned}[c]
  r & = 8 \\
  s & = 8 \\
  (\term{+} \enspace (\term{-} \enspace h \enspace 8) \enspace 1) & = 1024 \\
  (\term{+} \enspace (\term{-} \enspace w \enspace 8) \enspace 1) & = 256
\end{aligned}
\end{equation*}

The arithmetic simplification rule then produces:
\begin{equation*}
\begin{aligned}[c]
  i & = \term{[} 4 \term{,} 8 \term{,} h \term{,} w \term{]} \\
  n & = 4 \\
  k & = 4 \\
  c & = 8
\end{aligned}
\quad
\begin{aligned}[c]
  r & = 8 \\
  s & = 8 \\
  (\term{-} \enspace h \enspace 8) & = 1023 \\
  (\term{-} \enspace w \enspace 8) & = 255
\end{aligned}
\end{equation*}

And again:
\begin{equation*}
\begin{aligned}[c]
  i & = \term{[} 4 \term{,} 8 \term{,} h \term{,} w \term{]} \\
  n & = 4 \\
  k & = 4 \\
  c & = 8
\end{aligned}
\quad
\begin{aligned}[c]
  r & = 8 \\
  s & = 8 \\
  h & = 1031 \\
  w & = 263
\end{aligned}
\end{equation*}

The basic rule for variable replacement then produces:
\begin{equation*}
\begin{aligned}[c]
  i & = \term{[} 4 \term{,} 8 \term{,} 1031 \term{,} 263 \term{]} \\
  n & = 4 \\
  k & = 4 \\
  c & = 8
\end{aligned}
\quad
\begin{aligned}[c]
  r & = 8 \\
  s & = 8 \\
  h & = 1031 \\
  w & = 263
\end{aligned}
\end{equation*}

And shape inference is complete. The shape of the input has been inferred as $\term{[}4\term{,}8\term{,}1031\term{,}263\term{]}$.

\subsection{Discovering Shapes of Inputs}

The shape inference algorithm allows the shapes of unknown inputs to be inferred from other
information in the program. For example, consider the following example. The type signature of
\lstinline|max| determines the shape of the second input, given the shape of the first input. This
demonstrates that shape information can be inferred not only ``forward'' from input shapes to output
shapes, but also ``in reverse'' back to the inputs.
\begin{lstlisting}
fn(a : f32, b) {
  max(a, b)
}
\end{lstlisting}

Here the shape of \lstinline|b| is inferred to also be a scalar \term{f32}, due to the other input
being a scalar \term{f32} and the type signature for \lstinline{max}. The return type of the
function is also inferred to be \term{f32}.

\subsection{Partial Shape Information}

Sometimes there is not enough information to determine concrete shapes for all parts of a
program. Our shape inference algorithms allows partial shape information to be discovered, for
example consider the following \lstinline|softmax| program:
\begin{lstlisting}
softmax = fn (x : f32 s) {
  y = reduce(max, -inf, x);
  z = exp(x - y);
  z / reduce(+, 0f, z)
}
\end{lstlisting}

No shape information is given in the type of the input (it is just denoted \lstinline|s|), however
type information can be determined from the operations contained within the program. Running shape
inference on this program yields the following type for the \lstinline|softmax| function:
\begin{equation*}
\term{f32} \enspace \term{[}a\term{]} \term{@} b  \enspace \rightarrow \enspace \term{f32} \enspace \term{[}a\term{]} \term{@} b
\end{equation*}

Shape inference has discovered, based on the type signatures of the operations within the program,
that it can accept any input with a rank of at least one (the first dimension is denoted $a$) and
arbitrary outer dimensions (denoted $b$). Furthermore, the output of the function is also inferred
to have the same shape as the input.

\section{Introduction}
\label{sec:introduction}

Deep Learning models can be viewed as constrained functional programs on tensor domains, which only
permit side effects or updates for certain types of models, and usually only during training.
Tensors have a type and a shape: they are rectangular domains of simple element types. Requiring the
programmer to deal with these shapes can add significant complexity to a language.

This paper describes shape and rank inference for tensors in a language we call Axon, which aims to
make shapes implicit and inferred, in a similar manner to how types are implicit and inferred in
many functional programming languages. Axon allows the individual dimensions of a tensor to be
expressions consisting of symbolic variables, constants, and arithmetic operators, allowing complex
constraints on shapes to be expressed. Furthermore, shapes can be expressed as either a sequence of
these dimension expressions, as a symbolic variable, or as an appending the dimensions of several
other shapes. This allows rank inference to be expressed using this shape appending.

An inference algorithm is also presented which infers information about shapes, from both shape
information provided by the programmer and the structure of the program, via a set of rules for
built-in operators. Our system allows the shapes involved in complex Deep Learning problems to be
automatically inferred without need for the programmer to express the shape constraints or give
concrete shapes which are potentially unknown until the shape of the inputs to the program are
known. Shape information can be inferred from known shapes, or partial shapes, of program inputs and
from the structure of the program itself. Automated error detection also aids programmers by
identifying shape mismatches between the allowed shapes to an operation and inferred shapes.

The rest of the paper is structured as follows. Section~\ref{sec:language} outlines the Axon
language. Section~\ref{sec:syntax} presents the syntax of shape expressions,
section~\ref{sec:constraint-generation} describes how standard Hindley-Milner type
inference~\cite{Hindley69,Milner78} is used to generate a set of shape constraints for the program
and section~\ref{sec:shape-solver} describes the algorithm used to solve sets of these shape
constraints. Section~\ref{sec:examples} provides examples of shape inference in action and
section~\ref{sec:example-graphs} provides larger examples for commonly seen Deep Learning
graphs. Section~\ref{sec:related-work} discusses related work and section~\ref{sec:conclusions}
presents our conclusions.

This paper makes the following contributions:
\begin{itemize}
\item
  A taxonomy of the kinds of shapes encountered in Deep Learning graphs, which we use to guide the
  design of the language.
\item
  A functional programming language, which we call Axon. It allows a programmer to specify symbolic
  shapes for input and output tensors for a graph, which permit arithmetic expressions for the
  dimensions, and allows for rank inference by expressing shapes as the composition of sub-shapes.
\item
  A constraint based solver for inferring shape information throughout a program, based on a set of
  shape inference rules. Initial constraints are generated based on the structure of a graph, and
  the rules are used to reduce these constraints to a fixed point, where either all shapes are known
  rank with constant dimensions, or are partially unknown allowing for runtime variable shape.
\end{itemize}

\section{The Axon Language}
\label{sec:language}

The goal of Axon is to provide a language for easily writing Deep Learning graphs, in a target
agnostic manner. It provides a high-level functional language for representing Deep Learning graphs,
through the use of a suite of built-in deep learning operators. The Axon compiler includes a type
inference algorithm, so that explicit types do not need to be specified by the
programmer. Furthermore, an inference algorithm for shapes also removes the need for explicit tensor
shapes to be specified by the programmer.

A simple graph that compute the pointwise sum and product of two inputs can be expressed in Axon as
follows:
\begin{lstlisting}
f = fn(x, y) {
  a = x + y;
  b = x * y;
  a, b
}
\end{lstlisting}

This program describes a graph called $f$ which takes two inputs $x$ and $y$ and produces two
outputs $a$ and $b$. In Axon, all variables are tensors. These have an element type (with the usual
primitive types including fixed size integers and floating point data types, and a shape.

Axon allows program decomposition into sub-graphs as in the following example:
\begin{lstlisting}
g = fn(a, b) {
  a + b
};
f = fn(a, b) {
  g(a, b) + a;
}
\end{lstlisting}
This allows a program to be decomposed into smaller modular parts and reused. The type and shape
inference algorithm allow for polymorphic types and shapes. For example, in the above graph, if $g$
were called from multiple places it could be called with different types and shapes.

Axon includes a suite of operators for performing Deep Learning computations. These include matrix
multiplication, convolution, point-wise arithmetic operators, loops and more. The language includes
a type specification for each built-in, which is used to infer types and shapes throughout the
graph. Examples of these are given in Section~\ref{sec:example-shape-constraints}.

\section{Introduction}

\section{Related Work}
\label{sec:related-work}

Roesch~et~al.~\cite{Roesch2018RelayAN} describe RelayIR -- an intermediate representation for deep
learning. Similarly to Axon, RelayIR is a functional statically typed language aimed at deep
learning problems. In contrast to our work, it uses a dependent type system to allow shapes to be
represented by expressions. Our approach instead extends the language of shapes to include
arithmetic expressions, shape appending and other constructs in order to allow complex shapes to be
represented.

As far as we are aware there are no other related works for deep learning in the literature, that
use a constraint based methodology to infer shape information in graph programs.

Frameworks, such as TensorFlow~\cite{Tensorflow2025Whitepaper} and Pytorch~\cite{PytorchNIPS2019}
can propagate shape information at graph construction time but do not use an approach that allows
shape information to be inferred from outputs to inputs.

The MAGICA project~\cite{Joisha2006MATLAB} explores shape inference for MATLAB programs, using a symbolic
evaluation approach, based on the algebraic properties of MATLAB. They use it to identify redundant
shape array checks in MATLAB programs in order to reduce execution overheads at runtime.

Garg~et~al.~\cite{GargARRAY19} describe an approach that uses runtime information to perform shape
inference on APL-like languages for just-in-time compilation. In contrast, our inference algorithm
can be applied both statically at compile time, to infer potentially partial shape information, and
again at runtime (if required) once the shapes of graph inputs are known.

\section{Shape Solver}
\label{sec:shape-solver}

Type inference, as described in Section~\ref{sec:constraint-generation} produces a set of shape
constraints of the form:
\begin{equation*}
  C = \{ \dot{s}_1 = \dot{s}'_1, \ldots, \ldots \dot{s}'_n = \dot{s}'_n \}
\end{equation*}

The shape solver attempts to simplify this constraint set, and deduce concrete values for shapes
where possible. As shapes can depend on runtime values, the shape solver may not be able to
statically determine the value for all shape expressions. The solver uses an algorithm similar to
unification.

The $\text{solve-shapes}$ algorithm applies the following rules to the set of constraints until a
fixed point is reached. Rules are matched on the syntax of shapes, where $a$ denote shape variables,
$c$ denote integer constants, $s$ denotes shapes, $d$ denotes dimensions. $\dot{s}$ denotes both
shapes and dimensions. Two shapes are equivalent ($\equiv$) if their syntax is identical. If
\text{FAIL} is reached, the shape solver fails (as an invalid constraint has been found).

$C$ is the constraint set to be reduced. $C[\dot{s}/a]$ denotes the constraint set with all $a$
replaced with $\dot{s}$. $\text{fsv}(\dot{s})$ denotes the set of shape and dimension variables that
are free in $\dot{s}$.

\subsection*{Basic Rules}
These rules remove constraints that are tautologies, perform replacement of free shape variables
that have been determined and emit an error when different constants are equated to each other in a
constraint.
\begin{align*}
  C \cup \{ \dot{s}_1 = \dot{s}_2 \}
  & \Rightarrow C \enspace \text{if} \enspace \dot{s}_1 \equiv \dot{s}_2 \\
  C \cup \{ a = \dot{s} \}
  & \Rightarrow C[\dot{s}/a] \cup \{ a = \dot{s} \} \\
  & \qquad \text{if} \enspace a \in \text{fsv}(C) \enspace \text{and} \enspace
    a \notin \text{fsv}(\dot{s}) \\
  C \cup \{ c_1 = c_2 \}
  & \Rightarrow \text{FAIL} \enspace \text{if} \enspace c_1 \not\equiv c_2
\end{align*}

\subsection*{Reordering}
These rules reorder constraints with a constant on one side so that it is on the right hand
side. These rules remove the need for many symmetric rules later.
\begin{align*}
  C \cup \{ c = \dot{s} \}
  & \Rightarrow C \cup \{ \dot{s} = c \} \\
  C \cup \{ c = a \}
  & \Rightarrow C \cup \{ a = c \} \\
  C \cup \{ \term{*} = a  \}
  & \Rightarrow C \cup \{ a = \term{*} \} \\
  C \cup \{ \term{(op} \enspace d_1 \enspace d_2 \term{)} = a \}
  & \Rightarrow C \cup \{ a = \term{(op} \enspace d_1 \enspace d_2 \term{)} \} \\
  C \cup \{ \term{[} d_1 \term{,} \ldots \term{,} d_n \term{]} = a \}
  & \Rightarrow C \cup \{ a = \term{[} d_1 \term{,} \ldots \term{,} d_n \term{]} \} \\
  C \cup \{ s_1 \enspace \term{@} \enspace \ldots \enspace \term{@} \enspace s_n = a \}
  & \Rightarrow C \cup \{ a = s_1 \enspace \term{@} \enspace \ldots \enspace \term{@} \enspace s_n \} \\

\end{align*}

\subsection*{Unpacking}
These rules equate the shape expressions making up shapes of known rank.
\begin{align*}
  & C \cup \{ \term{[} d_1, \ldots \term{,} d_n \term{]} = \term{[} d'_1, \ldots \term{,} d'_n\term{]}
  \} \\
  & \Rightarrow C \cup \{ d_1 = d'_1, \ldots, d_n = d'_n \} \\
  & C \cup \{ \term{[} d_1, \ldots \term{,} d_n \term{]} = \term{[} d'_1, \ldots \term{,} d'_m
    \term{]} \} \\
  & \Rightarrow \text{FAIL} \enspace \text{if} \enspace n \neq m
\end{align*}

\subsection*{Empty shapes}
These rules assign the empty shape to shapes appended either side of a shape of known rank, where
possible.
\begin{align*}
  & C \cup \{ \term{[} d_1, \ldots \term{,} d_n \term{]} = s_1 \enspace \term{@} \ldots
  \term{@} \enspace s_j \enspace \term{@} \enspace \term{[} d'_1, \ldots \term{,} d'_n \term{]} \\
  & \enspace\enspace\enspace \term{@} \enspace s_{j+1} \enspace \term{@} \ldots \term{@} \enspace s_m \} \\
  & \Rightarrow C \cup \{ \term{[} d_1, \ldots \term{,} d_n \term{]} = \term{[} d'_1, \ldots \term{,}
  d'_n \term{]}, s_1 = \term{[]}, \ldots, s_m = \term{[]} \} \\
  & C \cup \{ s_1 \enspace \term{@} \ldots \term{@} \enspace s_j \enspace \term{@} \enspace
  \term{[} d'_1, \ldots \term{,} d'_n \term{]} \enspace \\
  & \enspace\enspace\enspace \term{@} \enspace s_{j+1} \enspace \term{@} \ldots \term{@} \enspace s_m = \term{[} d_1, \ldots \term{,} d_n \term{]} \} \\
  & \Rightarrow C \cup \{ \term{[} d_1, \ldots \term{,} d_n \term{]} = \term{[} d'_1, \ldots \term{,} d'_n \term{]},
  s_1 = \term{[]}, \ldots, s_m = \term{[]} \}
\end{align*}

\subsection*{Append Unpacking}
These rules equate dimension expressions that are at the start and end of appended shapes.
\begin{align*}
  & C \cup \{ s \enspace \term{@} \enspace \term{[} d_1 \term{,} \ldots d_n \term{]} =
    s' \enspace \term{@} \enspace \term{[} d'_1 \term{,} \ldots d'_m \term{]} \} \\
  & \Rightarrow C \cup \{ s \enspace \term{@} \enspace \term{[} d_1 \term{,} \ldots d_{n-1} \term{]} =
    s' \enspace \term{@} \enspace \term{[} d'_1 \term{,} \ldots d'_{m-1} \term{]}, d_n = d'_m \} \\
 & C \cup \{ \term{[} d_1 \term{,} \ldots d_n \term{]} \enspace \term{@} \enspace s =
    \term{[} d'_1 \term{,} \ldots d'_m \term{]} \enspace \term{@} \enspace s' \} \\
  & \Rightarrow C \cup \{ \term{[} d_2 \term{,} \ldots d_n \term{]} \enspace \term{@} \enspace s =
    \term{[} d'_2 \term{,} \ldots d'_m \term{]} \enspace \term{@} \enspace s', d_1 = d'_1  \}
\end{align*}

\subsection*{Append Dimensionality}
These rules check for shape appends whose minimum rank is greater than a concrete shape that they
are equal to.
\begin{align*}
  & C \cup \{ s_1 \enspace \term{@} \enspace \ldots \enspace \term{@} \enspace s_n =
    \term{[} d_1 \term{,} \ldots \term{,} d_m \term{]} \} \\
  & \Rightarrow \text{FAIL} \enspace \text{if number of dimensions in $s_1, \ldots, s_n > m$ } \\
  & C \cup \{ \term{[} d_1 \term{,} \ldots \term{,} d_m \term{]} =
    s_1 \enspace \term{@} \enspace \ldots \enspace \term{@} \enspace s_n \} \\
  & \Rightarrow \text{FAIL} \enspace \text{if number of dimensions in $s_1, \ldots, s_n > m$ }
\end{align*}

\subsection*{Simplification}
These rules use the $\text{simp}$ function to simplify shapes in the constraints.
\begin{align*}
  & C \cup \{ \dot{s}_1 = \dot{s}_2 \} \\
  & \Rightarrow C \cup \{ \dot{s}'_1 = \dot{s}'_2 \} \\
  & \qquad \text{where} \enspace
      \dot{s}'_1 = \text{simp}(\dot{s}_1), \enspace
      \dot{s}'_2 = \text{simp}(\dot{s}_2) \enspace \\
  & \qquad \text{if} \enspace
      \dot{s}_1 \neq \dot{s}_1' \enspace \text{or} \enspace \dot{s}_2 \neq \dot{s}'_2 \\
\end{align*}

The $\text{simp}$ function takes a shape, or a shape expression, and returns its canonical
form. This algorithm is recursive over the syntax of shapes and shape expressions. It simplifies
shapes, by computing operators of any constants, ordering variables by name inside arithmetic
operators, removing empty shapes in shape appends, and computing shape appends where possible.

\subsection*{Arithmetic Simplification}
These rules simplify arithmetic expressions.
\begin{align*}
  C \cup \{ \term{(*} \enspace d \enspace c_1 \term{)} = c_2 \}
  & \Rightarrow C \cup \{ s = \frac{c_2}{c_1} \} \enspace \text{if} \enspace c_2 \bmod c_1 = 0 \\
  C \cup \{ \term{(+} \enspace d \enspace c_1 \term{)} = c_2 \}
  & \Rightarrow C \cup \{ s = c_2 - c_1 \} \enspace \text{if} \enspace c_2 - c_1 > 0 \\
  C \cup \{ \term{(-} \enspace d \enspace c_1 \term{)} = c_2 \}
  & \Rightarrow C \cup \{ s = c_1 + c_2 \}
\end{align*}

\subsection*{Partial Expression Simplification}
This rule allows unsolvable expressions that equate to a constant to be replaced within other
constraints.
\begin{align*}
  & C + \{ (\term{op} \enspace s_1 \enspace s_2) = c \} \\
  & \Rightarrow C[c/(\term{op} \enspace s_1 \enspace s_2)] +
  \{ (\term{op} \enspace s_1 \enspace s_2) = c \} \\
  & \qquad \text{if} \enspace (\term{op} \enspace s1 \enspace s2) \in C
\end{align*}

\section{Tensor Shapes}
\label{sec:syntax}

In this section we describe the shapes of tensors in the Axon
language. Section~\ref{sec:taxonomy-of-shapes} outlines a general taxonomy of shapes we use to guide
our design and discussion. These are a general overview of the sorts of shapes that appear in Deep
Learning graphs. In Section~\ref{sec:syntax-of-shapes} we give the concrete syntax of shapes in the
Axon. Examples of shape constraints are given in Section~\ref{sec:example-shape-constraints}, which
show how shape constraints are generated from the syntactic structure of a graph. The algorithm for
then solving these constraints is given later, in Section~\ref{sec:shape-solver}

\subsection{Taxonomy of Shapes}
\label{sec:taxonomy-of-shapes}

Here we specify the general taxonomy of shapes we use to guide the design of our language. This
taxonomy categorizes the possible tensor shapes based on their properties and constraints.

\begin{itemize}
\item
  Concrete shapes

  These are shapes of statically known constant rank whose dimensions are all statically known
  constants.

  For example \lstinline|[]| denotes a scalar, \lstinline|[256]| denotes a one dimensional tensor of
  256 elements, and \lstinline|[1024, 256]| denotes a two dimension tensor.

\item
  Symbolic dimensions

  These extend concrete shapes by allowing individual dimensions to be symbolic expressions. The
  rank of the shape is still statically known.

  For example \lstinline|[x,y,16]| is a three dimensional tensor whose inner two dimensions are
  symbolic variables. At runtime, \lstinline|x| and \lstinline|y| will be known, but at compile time
  their value is unknown.

  \lstinline|[x*2,y,16]| is the shape that results from the concatenation of two tensors with the
  previous shape. Here the outermost dimension is a symbolic expression.
\item
  Symbolic

  These extend symbolic dimension shapes by allowing the for unknown rank.

  For example \lstinline|s @ [x,y]| denotes a tensor with at least two dimensions. \lstinline|s| is
  an arbitrary shape, and the \lstinline|@| operator is used to append \lstinline|s| to the shape
  \lstinline|[x,y]|.
\item
  Runtime dimensions

  Here the rank is static, but the individual dimensions can depend on the computation.
\item Runtime

  These extend runtime dimensions shapes by allowing the rank to also depend on computation.

\item
  Variable

  Dimension values differ across the tensor. This allows, for example, ragged batches to be
  represented.
\end{itemize}

\subsection{Syntax of Shapes}
\label{sec:syntax-of-shapes}

Figure~\ref{fig:shape-syntax} gives the syntax of shapes in Axon. The syntax of shapes is too
complex for a succinct grammar, therefore the grammar given in Figure~\ref{fig:shape-syntax}
generates some invalid shapes, which are further constrained by the shape construction rules given
in Figure~\ref{fig:shape-construction-rules}.

Shapes can be fixed, or variable, according to the dimension expressions they contain. The syntax of
these expressions includes positive integer constants, symbolic variables and arithmetic operations
and a ``variable dimension'' denoted \term{*}. The \term{*} dimension expression (with no arguments)
denotes a variable dimension, i.e. the number of elements in this dimension varies across the other
dimensions in the tensor. This allows, for example, batches of variable length sequences to be
represented. The shape of a batch of 10 variable length sequences of 3-dimensional vectors would be
denoted $\term{[10, *, 3]}$.

The shape append construct \term{@} allows for shapes to be expressed as many sub-shapes that are
appended to one another. This allows shape constraints such as $s \term{@} \term{[} n \term{]}$
which denotes a shape with a rank of at least one, where the outermost dimension has size
\term{n}. This is how Axon provides rank inference. We have found from our experience that this
append construct is able to represent the kinds of variable rank that we see in real Deep Learning
graphs.

\begin{figure}
  \begin{syntax}
    \nterm{Id}
       &=& [\nterm{Alpha} \enspace \term{\_}]
         [\nterm{Alpha} \enspace \nterm{Num} \enspace \term{\_}]^*
         & \text{Identifiers} \\
  \nterm{Alpha}
       &=& [\term{A}-\term{Z} \term{a}-\term{z}] \\
  \nterm{Num}
       &=& [\term{0}-\term{9}] \\
  \nterm{Shape}
       &=& \term{[]} & \text{Empty shape (scalar)} \\
    &\mid& \term{[} \nterm{Dimension} \big(\term{,} \enspace \nterm{Dimension}\big{)}^* \term{]}
         & \text{Shape with dimensions}\\
    &\mid& \nterm{Shape} \enspace \big( \term{@} \enspace \nterm{Shape} \big)^+
         & \text{Shape appending} \\
    &\mid& \nterm{Id} & \text{Shape variable} \\
  \nterm{Dimension}
       &=& \nterm{Num}^+ & \text{Dimension constant} \\
    &\mid& \term{*} & \text{Variable dimension} \\
    &\mid& \nterm{Id} & \text{Dimension variable} \\
    &\mid& \term{(} \term{+} \enspace \nterm{Dimension} \enspace \nterm{Dimension} \term{)}
         & \text{Addition} \\
    &\mid& \term{(} \term{-} \enspace \nterm{Dimension} \enspace \nterm{Dimension} \term{)}
         & \text{Subtraction} \\
    &\mid& \term{(} \term{*} \enspace \nterm{Dimension} \enspace \nterm{Dimension} \term{)}
         & \text{Multiplication} \\
    &\mid& \term{(} \term{/} \enspace \nterm{Dimension} \enspace \nterm{Dimension} \term{)}
         & \text{Division (floored)}
  \end{syntax}
  \caption{The syntax of shapes}
  \label{fig:shape-syntax}
\end{figure}

\begin{figure}
  \begin{typerules}
  \typerule{shape}
    {\forall i \in [1, n] \big( \vdash s_i : \nterm{Dimension} \big)}
    {\vdash \term{[} s_1 \term{,} \ldots \term{,} s_n \term{]} : \term{Shape}} \\
  \typerule{shape-append}
    {\forall i \in [1, n] \big( \vdash s_i : \nterm{Shape} \big)}
    {\vdash s_1 \enspace \term{@} \enspace \ldots \enspace \term{@} \enspace s_n : \term{Shape}} \\
  \typerule{shape-id}
    {s \in \nterm{Id}}
    {\vdash s : \term{Shape}} \\
  \typerule{dimension-constant}
    {d \in \nterm{Num}^+}
    {\vdash d : \term{Dimension}} \\
  \typerule{dimension-variable}
    { }
    {\vdash \term{*} : \term{Dimension}} \\
  \typerule{dimension-id}
    {d \in \nterm{Id}}
    {\vdash d : \term{Dimension}} \\
  \typerule{dimension-op}
    {\vdash d_1 : \term{Dimension} \and \vdash d_2 : \term{Dimension} \\ \and op \in \{+, -, *, /\}}
    {\vdash \term{(} op \enspace d_1 \enspace d_2 \term{)} : \term{Dimension}} \\
  \end{typerules}
  \caption{The rules for constructing shapes}
  \label{fig:shape-construction-rules}
\end{figure}

\subsection{Example Shape Constraints}
\label{sec:example-shape-constraints}

The following gives the types of commonly encountered deep learning operators. These type signatures
are used by the shape inference to infer information about the shapes of inputs and outputs and to
check that inputs and outputs to these operators are of a valid shape. When a graph is constructed,
rules such as these are used to create a set of shape constraints that can then be solved to infer
shape information about a graph. More details of how shape constraints are generated is given in
Section~\ref{sec:constraint-generation}, and the algorithm used to solve them is given in
Section~\ref{sec:shape-solver}.

\begin{itemize}
\item Concatenation \\
  $(\tau [d] @ s, \tau [d] @ s) \rightarrow \tau [(* \enspace d \enspace 2)] @ s$ \\
  This function takes two tensors which have shape $[d_1] @ s$ and $[d_2] @ s$, and concatenates
  them along their outermost dimension. This results in an output of shape $[(+ \enspace
    d_1 \enspace d_2)] @ s$. The outermost dimension of the output has a size which is equal to the
  sum of the outermost dimensions of the inputs, and the rest of the dimensions are arbitrary but
  match the inputs.
\item Matrix multiplication \\
  $(\tau s @ [d_1, d_2], \tau s @ [d_2, d_3]) \rightarrow \tau s @ [d_1, d_3]$ \\
  This function takes a tensor of shape $s @ [d_1, d_2]$ and a tensor of shape $s @ [d_1, d_2]$ and
  produces an output tensor of shape $s @ [d_1, d_3]$. $s$ is the shape of the batch (i.e. the
  matrix multiplication operates on the outer two dimensions of the inputs, and the outer dimensions
  must be equal, and are denoted $s$). The outer dimensions of the output shape are $d_1, d_3$,
  which expresses the fact that the result has the shape of the first dimension of the first input,
  and the second dimension of the second input. Finally, the inputs share dimension expression $d_2$
  which expresses the need for the outer and inner dimensions of the two matrices to be equal.
\item 2D Convolution in NCHW format \\
  $(\tau [n,c,h,w], \tau [k,c,r,s]) \rightarrow \tau [n, c, 1 + h - r, 1 + w - s]$ \\
  The shape of the output of a convolution is an arithmetic function of the shape of the inputs.
\end{itemize}