\section{Introduction}
\def\mathfrak{A}{\mathfrak{A}}
This work is inspired by the classical commutation relations in quantum physics such as the canonical commutation relation for bosonic systems \cite{Gaa54b}, the canonical anticommutation relation for fermionic systems \cite{Gaa54a}, the deformed commutation relation for $q$-oscillators \cite{Bie89,Hel00,Hel05,Mac89}, and also the related commutation relations for rotation algebras \cite{Eff67,Rie81,Web13}. The study of these commutation relations has led to the development of elegant results under Algebra and Functional Analysis. See for instance, \cite[pp. 217--220]{Bra97}. An abstract setting for all commutation relations mentioned is to view them as specific manifestations of a defining relation of a unital associative algebra $\mathfrak{A}$ generated by two elements $A$, $B$. We can express this defining relation for $\mathfrak{A}$ in the form $AB=\sigma(BA)$ for some function $\sigma$ defined on the free unital associative algebra on $A,B$ into itself, and this function $\sigma$, the deformation or twisting map, can be considered as ``linear'' in the sense of a polynomial of degree one. We call the relation $AB=\sigma(BA)$ as a \emph{deformed} or  \emph{twisted commutation relation}.

A natural Lie algebra structure is induced in $\mathfrak{A}$ by the operation $\left[ P,Q\right] :=PQ-QP$ defined for any $P,Q\in\mathfrak{A}$. We are interested in the Lie subalgebra of $\mathfrak{A}$ generated by the same generators $A,B$, that is, the set of all Lie polynomials in $A,B$. Our goal is to obtain a characterization for any nonzero $P\in\mathfrak{A}$ to be a Lie polynomial in $A,B$, or alternatively, to determine a linear complement of the Lie subalgebra of $\mathfrak{A}$ generated by $A,B$. We shall call the problem of finding such a characterization (or the corresponding direct sum decomposition) as the \emph{Lie polynomial characterization problem} for the algebra $\mathfrak{A}$ (with respect to the presentation in which $A,B$ are generators). As we shall see, nontrivial forms of $AB=\sigma(BA)$ are inexpressible in terms of Lie algebra operations only, which means that the results from the theory of bases of finitely-generated and finitely-presented Lie algebras is not directly applicable, and so other  methods have to be devised in order to solve the Lie polynomial characterization problem.

The solution to the Lie polynomial characterization problem for a free algebra has resulted into an elegant theory that involves Hopf algebras. See for instance \cite[Theorem 1.4]{Reu93} in which the maps used in the Lie polynomial characterizations are precursors to establishing a Hopf algebra structure in the same algebra \cite[Proposition 1.10]{Reu93}. For algebras that are not free, we are interested in algebras with defining relations that are not (all) Lie polynomials. The first appearance of such a Lie polynomial characterization problem is \cite[Problem 12.14]{Ter11}, which is about the universal Askey-Wilson algebra $\Delta$. Initial progress was made in \cite{Can15}. However, the Lie polynomial characterization problem for $\Delta$ remains open. For $q$-deformed Heisenberg algebras with $q$ not a root of unity, the Lie polynomial characterization problem is solved \cite{Can17}, and this current paper serves as an initial step in generalizing \cite{Can17}. Another recent result is a solution to the Lie polynomial characterization problem for a toral algebra (an algebra related to a torus) in which the Fairlie-Odesskii algebra is embedded \cite{Can18}. The Lie polynomial characterization in the said toral algebra was used to prove that the Casimir element of the Fairlie-Odesskii algebra is not a Lie polynomial if the deformation parameter is not a root of unity \cite[Theorem 5.3]{Can18}.

The basis theorem for the algebra in this paper, Lemma~\ref{basisLem}, is a generalization of \cite[Corollary 4.5]{Hel05}, which is for the specific case of a $q$-deformed Heisenberg algebra. The proof of \cite[Corollary 4.5]{Hel05} makes use of gradations, degree functions, chain conditions on ideals, and explicit reordering formulae that involve $q$-special combinatorics \cite[Appendix C]{Hel00}. Our proof for the generalization, Lemma~\ref{basisLem}, is elementary and makes use only of \emph{the Diamond Lemma for Ring Theory} \cite[Theorem 1.2]{Ber78}. Our main result, Theorem~\ref{mainThm}, is a generalization of \cite[Theorem 5.8]{Can17}, but our proof of the former is also elementary compared to the proof of the latter from \cite{Can17}, which relies on the theory of the \emph{Lyndon-Shirshov basis} of a free Lie algebra, which can be found in sources such as the seminal works \cite{Shi53,Shi58}, and the excellent exposition \cite[Section 2.8]{Ufn95}.

\section{Preliminaries}
Let ${\mathbb F}$ denote a fixed but arbitrary field. Throughout, an algebra is assumed to be unital and associative. Given an algebra $\mathcal{A}$ over ${\mathbb F}$ with unity element $I_\mathcal{A}$, we use the convention that $U^0=I_\mathcal{A}$ for any $U\in\mathcal{A}$, and we denote $I_\mathcal{A}$ simply by $I$ if no confusion shall arise. Any subalgebra is assumed to contain the unity element. Thus, if a subalgebra $\mathcal{S}$ of a nontrivial algebra $\mathcal{A}$ is generated by only one element $X\in\mathcal{A}\backslash\{0,I\}$, we simply write $\mathcal{S}$ as the polynomial algebra ${\mathbb F}[X]\subseteq\mathcal{A}$. We also note that every algebra $\mathcal{A}$ has a Lie algebra structure induced by $\left[ U,V\right]:=UV-VU$ for all $U,V\in\mathcal{A}$. Throughout, whenever we refer to the Lie algebra structure on an algebra $\mathcal{A}$, we shall always mean that which is induced by the Lie bracket operation $\left[\cdot,\cdot\right]$ just mentioned. Let $X_1,X_2,\ldots,X_n\in\mathcal{A}$. If $\mathcal{K}$ is the Lie subalgebra of $\mathcal{A}$ generated by $X_1,X_2,\ldots,X_n$, then we call the elements of $\mathcal{K}$ the \emph{Lie polynomials in $X_1,X_2,\ldots,X_n$}. Given $U\in\mathcal{A}$, the linear map ${\rm ad}\  U:\mathcal{A}\rightarrow\mathcal{A}$ is defined by the rule $V\mapsto\left[ U,V\right]$. Also, for any linear map $\varphi$ whose domain and codomain are equal, and for any nonnegative integer $n$, by $\varphi^n$ we mean composition of $\varphi$ with itself $n$ times, where $\varphi^0$ is interpreted as the identity linear map $\texttt{id}$. We interpret the composition of an empty collection of linear maps also as $\texttt{id}$.

Given a nonempty finite set $\mathcal{S}$, we denote by $\left<\mathcal{S}\right>$ the free monoid on $\mathcal{S}$. We call the elements of $\left<\mathcal{S}\right>$ the \emph{words} on $\mathcal{S}$. The multiplicative operation on $\left<\mathcal{S}\right>$ is the concatenation product, and we call the multiplicative identity in $\left<\mathcal{S}\right>$ the \emph{empty word}. By ${\mathbb F}\left<\mathcal{S}\right>$ we mean the free algebra generated by $\mathcal{S}$. The set $\left<\mathcal{S}\right>$ is a basis for the vector space ${\mathbb F}\left<\mathcal{S}\right>$, and the empty word in $\left<\mathcal{S}\right>$ is the unity element of the algebra ${\mathbb F}\left<\mathcal{S}\right>$.

Consider a two-element generating set $\mathcal{X}:=\{A,B\}$. Given a function  $\sigma:{\mathbb F}\left<\genset\right>\rightarrow{\mathbb F}\left<\genset\right>$, denote by $\mathcal{I}_\sigma$ the (two-sided) ideal of ${\mathbb F}\left<\genset\right>$ generated by $AB-\sigma(BA)$. We call the defining relation $AB=\sigma(BA)$ of the quotient algebra ${\mathbb F}\left<\genset\right>/\mathcal{I}_\sigma$ as \emph{the $\sigma$-twisted commutation relation (satisfied by $A,B$)}. If there exist $m,b\in{\mathbb F}$ with $m\neq 0$ such that $\sigma(X)=mX+bI$, we say that the commutation relation $AB=\sigma(BA)$ is \emph{linearly twisted (with parameters $m,b$)}. We may also refer to $m$ as the \emph{slope parameter} of the twisting map $\sigma$. We further exclude the case $m=1$, where ${\mathbb F}\left<\genset\right>/\mathcal{I}_\sigma$ becomes the \emph{Heisenberg-Weyl algebra} \cite[Definition 2.1]{Goo15} over ${\mathbb F}$. Algebraic aspects such as reordering and normal forms in the Heisenberg-Weyl algebra are excellently discussed in \cite[Section 1.2]{Hel02} and also in \cite{Bla03}. However, the Lie polynomial characterization problem in the Heisenberg-Weyl algebra has a trivial solution. The set of all Lie polynomials in $A,B$ is the three-dimensional Lie algebra with basis $A,B,\lbrack A,B\rbrack$, and the corresponding commutator table for this basis is completely determined by the relations $\left[ A,\lbrack A,B\rbrack\right]=0=\left[\lbrack A,B\rbrack, B\right]$.

From this point onward, assume that $AB=\sigma(BA)$ is linearly twisted with parameters $m,b$. By a simple result from the theory of free Lie algebras \cite[Lemma 1.7]{Reu93}, the generator $AB-\sigma(BA)$ of $\mathcal{I}_\sigma$ is not an element of the free Lie algebra on $\mathcal{X}$. That is, the defining relation of ${\mathbb F}\left<\genset\right>/\mathcal{I}_\sigma$ cannot be expressed in terms of Lie algebra operations only. Thus, the existing elegant and powerful machinery for computing a basis for finitely-presented Lie algebras, such as those in \cite{Ger96,Ger97a,Ger97b}, are not directly applicable, but in this work we present an elementary approach to the Lie polynomial characterization problem in  ${\mathbb F}\left<\genset\right>/\mathcal{I}_\sigma$.

\section{A reduction system for ${\mathbb F}\left<\genset\right>/\mathcal{I}_\sigma$ in three generators}

 As will become apparent in the succeeding parts of this section, we have good reason to consider properties of ${\mathbb F}\left<\genset\right>/\mathcal{I}_\sigma$ that are related to an algebra with three generators. Let $\mathcal{Y}=\{A,B,C\}$, and denote by $\mathcal{J}_\sigma$ the  ideal of ${\mathbb F}\left<\newgenset\right>$ generated by $AB-mBA-bI$ and $C-AB+BA$. An immediate consequence is that there exists an algebra isomorphism ${\mathbb F}\left<\genset\right>/\mathcal{I}_\sigma\rightarrow{\mathbb F}\left<\newgenset\right>/\mathcal{J}_\sigma$ such that $A\mapsto A$, $B\mapsto B$, and $\lbrack A,B\rbrack\mapsto C$.

 Our goal in this section is to determine a basis for $\F\newmonoid/\newidcomm_\sigma$ using the Diamond Lemma for Ring Theory \cite[Theorem 1.2]{Ber78}. We establish a sufficient condition for some elements of $\F\newmonoid/\newidcomm_\sigma$ to form a basis for $\F\newmonoid/\newidcomm_\sigma$, by invoking the Diamond Lemma. To establish this sufficient condition, we recall terminology and notation from \cite[Section 1]{Ber78}. We introduce each of such notions as the need for it arises. As a first step, we give an alternative presentation of $\F\newmonoid/\newidcomm_\sigma$ that will turn out to be more suited to the Diamond Lemma.

\begin{proposition}\label{presProp} The algebra ${\mathbb F}\left<\newgenset\right>/\mathcal{J}_\sigma$ has a presentation with generators $A,B,C$ and relations
\begin{eqnarray}
AB & = &\frac{mC-bI}{m-1},\label{redrel1}\\
AC & = & mCA,\label{redrel2}\\
BA & = &\frac{C-bI}{m-1},\label{redrel3} \\
CB & = &  mBC,\label{redrel4}\\
BC^kA & = & \frac{C^{k+1}-bC^k}{m^k(m-1)},\quad\quad k\in\{1,2,\ldots\}.\label{redrel5}
\end{eqnarray}
\end{proposition}
\def\zeta{\zeta}
\def\xi{\xi}
\begin{proof} View the left-hand side $AB$ and the right-hand side $\frac{mC-bI}{m-1}$ of \eqref{redrel1} as elements of the free algebra ${\mathbb F}\left<\newgenset\right>$, and define $\xi_1\in{\mathbb F}\left<\newgenset\right>$ as the difference $AB-\frac{mC-bI}{m-1}$. In a similar manner, define $\xi_2,\  \xi_3,\  \xi_4\in{\mathbb F}\left<\newgenset\right>$ for the relations \eqref{redrel2} to \eqref{redrel4}, respectively. As for the countably infinite relations \eqref{redrel5}, for each positive integer $k$, we define $\xi_5(k)\in{\mathbb F}\left<\newgenset\right>$ to be the left-hand side of \eqref{redrel5} minus the right-hand side. Also, denote the generators of $\mathcal{J}_\sigma$ by $\zeta_1:=AB-mBA-bI$ and $\zeta_2:=C-AB+BA$. Let $\mathcal{K}$ be the ideal of ${\mathbb F}\left<\newgenset\right>$ generated by
\begin{eqnarray}
\{\xi_1,\xi_2,\xi_3,\xi_4\}\cup\{\xi_5(k)\  :\  k\in\{1,2,\ldots\}\}.
\end{eqnarray}
We claim that $\mathcal{J}_\sigma=\mathcal{K}$. By routine computations, it can be shown that the relations
\begin{eqnarray}
\zeta_1 & = & \xi_1-m\xi_2,\label{eqId1}\\
\zeta_2 & = &-\xi_1+\xi_2,\label{eqId2}\\
\xi_1 & = &\frac{\zeta_1+m\zeta_2}{1-m},\label{eqId3}\\
\xi_2 & = &\left[ A,\zeta_1\right]+A\zeta_1-m\zeta_2A,\label{eqId4}\\
\xi_3 & = &\frac{\zeta_1+\zeta_2}{1-m},\label{eqId5}\\
\xi_4 & = &\left[ \zeta_1,B\right]+\zeta_2B-mB\zeta_2,\label{eqId6}\\
\xi_5(k) & = & m^{-k}\xi_3-\sum_{i=1}^km^{i-1-k}C^{k-i}\xi_4C^{i-1}A,\   (k\in\{1,2,\ldots\}),\label{eqId7}
\end{eqnarray}
hold in ${\mathbb F}\left<\newgenset\right>$. We see from \eqref{eqId1} and \eqref{eqId2} that the generators $\zeta_1$, $\zeta_2$ of $\mathcal{J}_\sigma$ are elements of $\mathcal{K}$, and so $\mathcal{J}_\sigma\subseteq\mathcal{K}$. By \eqref{eqId3} to \eqref{eqId6}, the generators $\xi_1,\xi_2,\xi_3,\xi_4$ of $\mathcal{K}$ are in $\mathcal{J}_\sigma$. Then by \eqref{eqId7}, all the other generators $\xi_5(k)$ for all $k$ are also elements of $\mathcal{J}_\sigma$. Then $\mathcal{K}$ is contained in $\mathcal{J}_\sigma$, and hence $\mathcal{K}=\mathcal{J}_\sigma$. Therefore, ${\mathbb F}\left<\newgenset\right>/\mathcal{J}_\sigma={\mathbb F}\left<\newgenset\right>/\mathcal{K}$.
\end{proof}
We use the presentation of ${\mathbb F}\left<\newgenset\right>/\mathcal{J}_\sigma$ given in Proposition~\ref{presProp} to construct a \emph{reduction system in ${\mathbb F}\left<\newgenset\right>$} (i.e., a subset of $\left<\newgenset\right>\times{\mathbb F}\left<\newgenset\right>$). Let
\begin{eqnarray}
\alpha & := & \left(AB,\frac{mC-bI}{m-1}\right),\\
\beta & := & \left(AC,mCA\right),\\
\gamma & := & \left(BA,\frac{C-bI}{m-1}\right),\\
\delta & := & \left(CB,mBC\right),\\
\varepsilon(k) & := & \left(BC^kA , \frac{C^{k+1}-bC^k}{m^k(m-1)}\right),\quad k\in\{1,2,\ldots\}.
\end{eqnarray}
Then $\Gamma:=\{\alpha,\beta,\gamma,\delta\}\cup\{\varepsilon(k)\  :\  k\in\{1,2,\ldots\}\}$ is a reduction system in ${\mathbb F}\left<\newgenset\right>$ for the presentation of ${\mathbb F}\left<\newgenset\right>/\mathcal{J}_\sigma$ given in Proposition~\ref{presProp}. We use Bergman's notation for elements of a reduction system, which is as follows. Given an element $\mu$ of a reduction system, we denote the word that serves as the first coordinate of $\mu$ as $W_\mu$ and the second coordinate, a linear combination of words, by $f_\mu$. In the reduction system $\Gamma$ for example, we have $W_{\varepsilon(k+1)}=BC^{k+1}A$ and $f_\beta=mCA$.

If $S$ is a reduction system, a word $W\in\left<\newgenset\right>$ is said to be \emph{$S$-irreducible} whenever $W_\mu$ is not a subword of $W$ for any $\mu\in S$. It is routine to show that a word $W\in\left<\newgenset\right>$ is $\Gamma$-irreducible if and only if
\begin{eqnarray}
W=C^kA^l,\mbox{ or }\     W=B^lC^k,\label{iffirred}
\end{eqnarray}
for some nonnegative integers $k,l$.

We now recall the notion of ambiguities. Suppose $S$ is a reduction system in ${\mathbb F}\left<\newgenset\right>$. An \emph{overlap ambiguity} is a $5$-tuple
\begin{eqnarray}
\Phi=(\mu,\nu,L,X,R)\in S^2\times\left<\newgenset\right>^3\nonumber
\end{eqnarray}
such that
\begin{eqnarray}
W_\mu = LX,\quad W_\nu = XR.\label{overlapdef}
\end{eqnarray}
If instead of \eqref{overlapdef}, we have $W\mu=X$ and $W_\nu = LXR$, then $\Phi$ is called an \emph{inclusion ambiguity}. The term \emph{ambiguity} means either an overlap ambiguity or an inclusion ambiguity.

It is routine to show that the reduction system $\Gamma$ has no inclusion ambiguities, and we describe the inclusion ambiguities in the following. All overlap ambiguities that do not involve an element of $\Gamma$ of the form $\varepsilon(k)$ are the following:
\begin{eqnarray}
\Phi_1 & := & \left(\alpha,\gamma,A,B,A\right),\nonumber\\
\Phi_2 & := & \left(\beta,\delta,A,C,B\right),\nonumber\\
\Phi_3 & := & \left(\gamma,\alpha,B,A,B\right),\nonumber\\
\Phi_4 & := & \left(\gamma,\beta,B,A,C\right),\nonumber\\
\Phi_5 & := & \left(\delta,\gamma,C,B,A\right),\nonumber
\end{eqnarray}
while all the overlap ambiguities that depend on an integer parameter are:
\begin{eqnarray}
\Phi_6(k) & := & \left(\alpha,\varepsilon(k),A,B,C^kA\right),\nonumber\\
\Phi_7(k) & := & \left(\delta,\varepsilon(k),C,B,C^kA\right),\nonumber\\
\Phi_8(k) & := & \left(\varepsilon(k),\alpha,BC^k,A,B\right),\nonumber\\
\Phi_9(k) & := & \left(\varepsilon(k),\beta,BC^k,A,C\right).\nonumber
\end{eqnarray}

If $S$ is a reduction system on ${\mathbb F}\left<\newgenset\right>$, and if $L,R\in\left<\newgenset\right>$ and $\mu\in S$, then by a \emph{reduction} we mean a linear map $r_{L\mu R}:{\mathbb F}\left<\genset\right>\rightarrow{\mathbb F}\left<\genset\right>$ that fixes every element of the basis $\left<\newgenset\right>$ of ${\mathbb F}\left<\newgenset\right>$ except $LW_\mu R$, for which the rule $LW_\mu R\mapsto Lf_\mu R$ applies. For simpler notation, we write $r_{I\mu R}$, $r_{L\muI}$, $r_{I\muI}$ as $r_{\mu R}$, $r_{L\mu}$, $r_{\mu}$, respectively. An overlap ambiguity $\Phi=(\mu,\nu,L,X,R)$ is \emph{resolvable} if there exist compositions $\lambda,\rho$ of reductions such that
\begin{eqnarray}
\lambda(f_\mu R)=\rho(Lf_\nu).\label{overlapresolve}
\end{eqnarray}
We have an analogous definition for the \emph{resolvability of an inclusion ambiguity} in which \eqref{overlapresolve} is replaced by $\lambda(Lf_\mu R)=\rho(f_\nu)$.

For any positive integer $k$, the compositions
\begin{eqnarray}
a_k & := & r_{C^{k-1}\beta}\circ r_{C^{k-2}\beta C}\circ r_{C^{k-2}\beta C^2}\circ\cdots\circ r_{C\beta C^{k-2}}\circ r_{\beta C^{k-1}},\nonumber\\
b_k & := & r_{C^{k-1}\delta}\circ r_{C^{k-2}\delta C}\circ r_{C^{k-2}\delta C^2}\circ\cdots\circ r_{C\delta C^{k-2}}\circ r_{\delta C^{k-1}},\nonumber
\end{eqnarray}
of reductions for the reduction system $\Gamma$ shall be of importance in the succeeding lemmas in this section. We note that if $k=1$, we retain the subscript in $b_k=b_1$ so as not to cause confusion with the parameter $b$ of the linearly twisted commutation relation $AB=\sigma(BA)$. It is routine to show that for any positive integers $h,k,l$, the linear maps $a_k$ and $b_k$ satisfy the properties
\begin{eqnarray}
a_k(AC^l) & = & \left\{
    \begin{array}{ll}
      m^kC^kA, & l=k, \\
      AC^l, & l\neq k,
    \end{array}\right.\label{akbkprop1}\\
b_k(C^lB) & = & \left\{
    \begin{array}{ll}
      m^kBC^k, & l=k, \\
      C^lB, & l\neq k.
    \end{array}\right.\label{akbkprop2}\\
a_k(C^hA) & = & C^hA,\label{akbkprop3}\\
b_k(BC^h) & = & BC^h.\label{akbkprop4}
\end{eqnarray}

\begin{lemma}\label{ambigLem} All ambiguities of $\Gamma$ are resolvable.
\end{lemma}
\begin{proof} Consider the overlap ambiguity $\Phi_1$. To show that $\Phi_1$ is resolvable, we exhibit compositions $\lambda_1$, $\rho_1$ such that
\begin{eqnarray}
\lambda_1(f_\alpha A)=\rho_1(Af_\gamma).\label{reductionEq}
\end{eqnarray} We claim that $\lambda_1$ is the empty composition $\texttt{id}$ and that $\rho_1=r_\beta$. Since $r_\beta$ is a linear map,
\begin{eqnarray}
r_\beta(Af_\gamma)=\frac{r_\beta(AC)-br_\beta(A)}{m-1}.\label{resolve}
\end{eqnarray}
The key step is to examine the words in the right-hand side of \eqref{resolve} and determine if one of them is related to the element of the reduction system $\Gamma$ that appears in the subscript of the reduction involved. We see that $AC=W_\beta$, and so the right-hand side of \eqref{resolve} is a linear combination of $r_\beta(W_\beta)$ and $r_\beta(A)$. But since $AC$ and $A$ are linearly independent in ${\mathbb F}\left<\newgenset\right>$, we have $A\neq CA = W_\beta$, which means that $r_\beta$ fixes $A$, and by the definition of a reduction, we have $r_\beta(AC)=r_\beta(W_\beta)=f_\beta$. By these observations, \eqref{resolve} becomes
\begin{eqnarray}
r_\beta(Af_\gamma)=\frac{f_\beta-bA}{m-1}=\frac{mCA-bA}{m-1}=f_\alpha A=\texttt{id} (f_\alpha A),\nonumber
\end{eqnarray}
which implies that we have determined the solution $\lambda_1=\texttt{id}$, $\rho_1=r_\beta$ to \eqref{reductionEq}, as claimed. The ambiguities $\Phi_2,\ldots,\Phi_5,\Phi_6(k),\ldots,\Phi_9(k)$ can be shown to be resolvable by similar routine computations and arguments. We have an analogue of the equation  \eqref{reductionEq} for each $i\in\{2,3,\ldots,9\}$. We summarize the solutions in the following table.
\begin{center}
\begin{tabular}{|l|c|c|}
\hline
$i$ & $\lambda_i$  & $\rho_i$  \\
 \hline
1 & $\texttt{id}$ & $r_\beta$ \\
 \hline
2 & $r_{C\alpha}$ & $r_{\alpha C}$ \\
 \hline
3 & $r_{\delta}$  & $\texttt{id}$ \\
 \hline
4 & $\texttt{id}$  & $r_{\varepsilon(1)}$ \\
 \hline
5 & $r_{\varepsilon(1)}$ & $\texttt{id}$ \\
 \hline
6 & $\texttt{id}$ & $a_k\circ a_{k+1}$ \\
 \hline
7 & $r_{\varepsilon(k+1)}$ & $\texttt{id}$ \\
 \hline
8 & $b_k\circ b_{k+1}$ & $\texttt{id}$ \\
 \hline
9 & $\texttt{id}$ & $r_{\varepsilon(k+1)}$ \\
 \hline
\end{tabular}
\end{center}
More precisely, for each $i\in\{1,2,3,\ldots,9\}$, write the overlap ambiguity as either  $\Phi_i=(\mu,\nu,L,X,R)$ if $i<6$, or $\Phi_i(k)=(\mu,\nu,L,X,R)$ if $i\geq 6$ in which case one of $\mu$, $\nu$ depends on $k$. Then the above table shows the compositions $\lambda_i$, $\rho_i$ of reductions such that $\lambda_i(f_\mu R)=\rho_i(Lf_\nu)$, and this can be verified to be true for each $i\geq 2$ by routine computations. We only note here that in the computations for the case $i\in\{6,7,8,9\}$, the relations \eqref{akbkprop1} to \eqref{akbkprop4} are of significance. By the routine computations just described, all overlap ambiguities of $\Gamma$ are resolvable, and since there is no inclusion ambiguity, the proof is now complete.
\end{proof}

\begin{lemma}\label{basisLem} The vectors
\begin{eqnarray}
C^kA^l,\    B^lC^k,\     C^k\quad\in\quad\F\newmonoid/\newidcomm_\sigma, & & (k\in\{0,1,2,\ldots\},\label{theBasis}\\
& & \  \   l\in\{1,2,\ldots\},)\nonumber
\end{eqnarray}
form a basis for $\F\newmonoid/\newidcomm_\sigma$.
\end{lemma}
\begin{proof} Let $\Psi:{\mathbb F}\left<\newgenset\right>\rightarrow\F\newmonoid/\newidcomm_\sigma$ be the canonical map. Observe that the images of all $\Gamma$-irreducible words under $\Psi$ are precisely the vectors \eqref{theBasis}. The only difference between the conditions on the exponents $k,l$ noticeable in comparing \eqref{iffirred} and \eqref{theBasis} is that we avoid the duplication $C^kA^l=I=B^lC^k$ when $l=0=k$. Enough has been established in this section for us to invoke the Diamond Lemma \cite[Theorem 1.2]{Ber78}. The only implication we need from \cite[Theorem 1.2]{Ber78} is that: if all the ambiguities of a reduction system $S$ are resolvable and if $\mathcal{K}$ is the ideal of ${\mathbb F}\left<\newgenset\right>$ generated by all  $W_\mu-f_\mu$ ($\mu\in S$), then the images of all the $S$-irreducible words under the canonical map ${\mathbb F}\left<\newgenset\right>\rightarrow{\mathbb F}\left<\newgenset\right>/\mathcal{K}$ form a basis for ${\mathbb F}\left<\newgenset\right>/\mathcal{K}$. Thus if $S=\Gamma$ and $\mathcal{K}=\mathcal{J}_\sigma$, then by Lemma~\ref{ambigLem}, the vectors \eqref{theBasis} form a basis for ${\mathbb F}\left<\newgenset\right>/\mathcal{J}_\sigma$.
\end{proof}

\section{The Lie subalgebra of $\F\monoid/\mathcal{I}_\sigma$ generated by $A,B$}
All computations in this section are in the quotient algebra $\F\monoid/\mathcal{I}_\sigma$. We use the basis \eqref{theBasis} of ${\mathbb F}\left<\newgenset\right>/\mathcal{J}_\sigma=\F\monoid/\mathcal{I}_\sigma$ in order to distinguish which elements of $\F\monoid/\mathcal{I}_\sigma$ are Lie polynomials in $A,B$, and which are not. First, we give generalized forms of the relations \eqref{redrel2} and \eqref{redrel4}.
\begin{proposition} For any positive integers $k,l$, the relations
\begin{eqnarray}
A^lC^k & = & m^{kl}C^kA^l,\label{reorderAC}\\
C^kB^l & = &  m^{kl}B^lC^k,\label{reorderBC}
\end{eqnarray}
hold in $\F\monoid/\mathcal{I}_\sigma$.
\end{proposition}
\begin{proof} Use induction on $k,l$.
\end{proof}

\begin{example}\label{tableEx1} The relations \eqref{reorderAC} and \eqref{reorderBC} are reordering formulae for two types of products in $\F\monoid/\mathcal{I}_\sigma$. One is a product of a power of $A$ and a power of $C$ (with positive exponents), and the second type is a product of a power of $C$ and a power of $B$ (also with positive exponents). The importance of these reordering formulae is that given some important expressions like those shown in the left-hand sides of the equations below, we can rewrite such expressions in terms of the basis \eqref{theBasis} of $\F\monoid/\mathcal{I}_\sigma$. That is, by using \eqref{reorderAC}, \eqref{reorderBC}, it is routine to show that the relations
\begin{eqnarray}
\left[ A,C^x\right] & = & (m^{x}-1)C^{x}A,\label{table1}\\
\left[ A,C^xA^y\right] & = & (m^{x}-1)C^{x}A^{y+1},\label{table2}\\
\left[ C^k,C^xA^y\right] & = & (1-m^{ky})C^{k+x}A^{y},\label{table3}\\
\left[ C^kA^l,C^xA^y\right] & = & (m^{lx}-m^{ky})C^{k+x}A^{l+y},\label{table4}\\
\left[ B,C^x\right] & = & (1-m^{x})BC^{x},\label{table5}\\
\left[ B,B^yC^x\right] & = & (1-m^{x})B^{y+1}C^{x},\label{table6}\\
\left[ C^k,B^yC^x\right] & = & (m^{xy}-1)B^{y}C^{k+x},\label{table7}\\
\left[ B^lC^k,B^yC^x\right] & = & (m^{ky}-m^{lx})B^{l+y}C^{k+x}.\label{table8}
\end{eqnarray}
hold in $\F\monoid/\mathcal{I}_\sigma$ for any positive integers $k,l,x,y$. The relations \eqref{table1} to \eqref{table8} shall be important in the proof of our main result.
\end{example}

\begin{example}\label{tableEx2} Reordering in $\F\monoid/\mathcal{I}_\sigma$ will turn out to be more complicated if all of $A,B,C$ appear in an expression that is to be reordered. More of the relations in the presentation of $\F\monoid/\mathcal{I}_\sigma$ in Proposition~\ref{presProp}, will have to be involved, and not just \eqref{redrel2} and \eqref{redrel4} which have been generalized into \eqref{reorderAC} and \eqref{reorderBC}. As an example, by routine computations involving \eqref{redrel1}, \eqref{redrel3}, \eqref{reorderAC}, \eqref{reorderBC}, it can be shown that the relations
\begin{eqnarray}
\left[ A, B^{y}C^{x}\right] & = & \frac{m^y-m^{-x}}{m-1}B^{y-1}C^{x+1} + \frac{m^{-x}-m}{m-1}bB^{y-1}C^{x} ,\label{tablei}\\
\left[ B, C^{x}A^{y}\right] & = & \frac{m^{-x}-m^y}{m-1}C^{x+1}A^{y-1} + \frac{1-m^{-x}}{m-1}bC^{x}A^{y-1},\label{tableii}
\end{eqnarray}
hold in $\F\monoid/\mathcal{I}_\sigma$ for any positive integers $x,y$. The relations \eqref{tablei} and \eqref{tableii} shall also be of significance in the proof of our main result.
\end{example}

The partial descriptions of reordering in $\F\monoid/\mathcal{I}_\sigma$ described in Examples~\ref{tableEx1} and \ref{tableEx2} are sufficient for our purpose of solving the Lie polynomial characterization problem for the algebra $\F\monoid/\mathcal{I}_\sigma$. We now proceed with some technical lemmas.

\begin{lemma}\label{equalexpLem} For any positive integer $n$, there exist $f,g\in\F[C]\subseteq\F\monoid/\mathcal{I}_\sigma$ such that
\begin{eqnarray}
(m-1)^nA^nB^n & = & (-1)^nb^nI + Cf,\label{divAB}\\
(m-1)^nB^nA^n & = & (-1)^nb^nI + Cg.\label{divBA}
\end{eqnarray}
\end{lemma}
\begin{proof} We claim that
\begin{eqnarray}
(m-1)^nA^nB^n & = & \prod_{i=1}^{n}(m^{i}C-bI),\label{equalAB}\\
(m-1)^nB^nA^n & = & \prod_{i=0}^{n-1}(m^{-i}C-bI).\label{equalBA}
\end{eqnarray}
We use induction on $n$. If $n=1$, then \eqref{equalAB}, \eqref{equalBA} are immediate consequences of \eqref{redrel1}, \eqref{redrel3}, respectively. Suppose that \eqref{equalAB}, \eqref{equalBA} hold for some positive integer $n$. Multiply both sides of \eqref{equalAB} by $A$ from the left and $B$ from the right to obtain
\begin{eqnarray}
(m-1)^nA^{n+1}B^{n+1} & = & \prod_{i=1}^{n}(m^{i}ACB-bAB).\label{inductAB}
\end{eqnarray}
Use \eqref{redrel2} to replace $AC$ in every occurence of $ACB$ in the right-hand side of \eqref{inductAB}. The result is $h\cdot AB$ where $h$ is a polynomial in $C$. We then use \eqref{redrel1} to replace $AB$ by a polynomial in $C$. After some adjustments in the scalar coefficients, we find that \eqref{equalAB} holds for $n+1$. By similar computations, \eqref{equalBA} can also be shown to hold for $n+1$, and this proves the claim. It is routine to show that the coefficient of $I$ in the right-hand side of \eqref{equalAB} is $(-1)^nb^n$. Thus, every other term in the right hand side of \eqref{equalAB} is a scalar multiple of a power of $C$ with exponent at least 1. By these observations, the right-hand side of \eqref{equalAB} is equal to $(-1)^nb^nI + Cf$ for some polynomial $f$ in $C$. This proves \eqref{divAB}. The relation \eqref{divBA} can be proven similarly.
\end{proof}

\begin{lemma}\label{tableLem} If $u,v,x,y$ are positive integers, then $\left[ C^uA^v,B^yC^x\right]$ is an element of the span of
\begin{eqnarray}
C^k,\  C^kA^l,\  B^lC^k,\quad\quad (k,l\in\{1,2,\ldots\}).\label{derLieBasis}
\end{eqnarray}
\end{lemma}
\begin{proof} We first consider the case $v\geq y$. The trick is that $\left[ C^uA^v,B^yC^x\right]$ can be written as
\begin{eqnarray}
\left[ C^uA^v,B^yC^x\right] = C^uA^{v-y}(A^yB^y)C^x-B^yC^{u+x}A^yA^{v-y},\label{vatleasty}
\end{eqnarray}
simply because $v-y$ is a nonnegative integer. Reorder the expressions $C^uA^{v-y}$ and $B^yC^{u+x}$ that appear in the right-hand side of \eqref{vatleasty} using \eqref{reorderAC},\eqref{reorderBC}. We obtain
\begin{eqnarray}
\left[ C^uA^v,B^yC^x\right] = m^{u(y-v)}A^{v-y}P-m^{-y(u+x)}QA^{v-y},\label{vatleasty2}
\end{eqnarray}
where $P=C^u(A^yB^y)C^x$ and $Q=C^{u+x}(B^yA^y)$. By Lemma~\ref{equalexpLem}, $P$ and $Q$ are both polynomials in $C$. Also, the conditions $u>0$ and $x>0$ imply that the powers of $C$ that appear in the linear combinations for $P$ and $Q$ as polynomials in $C$ have exponents that are each at least $1$. This further implies that the first term $m^{u(y-v)}A^{v-y}P$, in the right-hand side of \eqref{vatleasty2}, is a linear combination of products of powers of $A$ and of $C$, which can be reordered using \eqref{reorderAC} such that  $m^{u(y-v)}A^{v-y}P$ can now be written as a linear combination of
\begin{eqnarray}
C^kA^{v-y},\quad\quad (k\in\{1,2,\ldots\}).\label{CAvectors}
\end{eqnarray}
By inspection, $-m^{-y(u+x)}QA^{v-y}$, the second term in the right-hand side of \eqref{vatleasty2}, is also a linear combination of \eqref{CAvectors}, and by \eqref{vatleasty2}, so is $\left[ C^uA^v,B^yC^x\right]$. The vectors \eqref{CAvectors} are among \eqref{derLieBasis}, and this completes the proof for the case $v\geq y$. By similar computations and arguments, it can be shown that for the case $v<y$, the expression $\left[ C^uA^v,B^yC^x\right]$ is a linear combination of
\begin{eqnarray}
B^{y-v}C^k,\quad\quad (k\in\{1,2,\ldots\}),
\end{eqnarray}
which are also among \eqref{derLieBasis}.
\end{proof}

\begin{lemma}\label{CpowerLem} For any positive integer $k$, the relation
\begin{eqnarray}
(1-m^{k+1})C^{k+1}=(1-m^k)bC^k-\frac{\left(\left({\rm ad}\  B\right)\circ\left({\rm ad}\  C\right)^k\right)(A)}{(1-m)^{k-1}m^{-k}}\label{powerofC}
\end{eqnarray}
holds in $\F\monoid/\mathcal{I}_\sigma$.
\end{lemma}
\begin{proof} Using the relation \eqref{redrel2}, it can be shown by induction that the relation
\begin{eqnarray}
({\rm ad}\  C)^k (A) = (1-m)^k C^kA\label{adCA}
\end{eqnarray}
holds in $\F\monoid/\mathcal{I}_\sigma$ for each positive integer $k$. Applying ${\rm ad}\  B$ on both sides of \eqref{adCA}, we get
\begin{eqnarray}
(({\rm ad}\  B)\circ({\rm ad}\  C)^k) (A) = (1-m)^k (BC^kA-C^kAB).\label{adBadCA}
\end{eqnarray}
Using \eqref{reorderBC}, the right-hand side of \eqref{adBadCA} becomes $$(1-m)^kC^k(m^{-k}BA-AB),$$
in which $BA$ and $AB$ can be replaced by polynomials in $C$ using \eqref{redrel1}, \eqref{redrel3}. In the resulting equation, we solve for $(1-m^{k+1})C^{k+1}$, and the result is \eqref{powerofC}.
\end{proof}

\begin{theorem}\label{mainThm} If $m$ is not a root of unity, then the span of the basis vectors
\begin{eqnarray}
I,\   A^n,\  B^n,\quad\quad (n\geq 2)\label{compLieBasis}
\end{eqnarray}
of $\F\monoid/\mathcal{I}_\sigma$ is a linear complement of the Lie subalgebra of $\F\monoid/\mathcal{I}_\sigma$ generated by $A,B$.
\end{theorem}
\def\mathcal{K}{\mathcal{K}}
\def\mathcal{L}{\mathcal{L}}
\begin{proof} Denote by $\mathcal{L}$ the Lie subalgebra of $\F\monoid/\mathcal{I}_\sigma$ generated by $A,B$, and let $\mathcal{K}$ be the span of all vectors in the basis \eqref{theBasis} of $\F\monoid/\mathcal{I}_\sigma$ that are not in \eqref{compLieBasis}. That is, the vectors
\begin{eqnarray}
A,\  B,\  C^k,\  C^kA^l,\  B^lC^k,\quad\quad (k,l\in\{1,2,\ldots\}),\label{LieBasis}
\end{eqnarray}
form a basis for $\mathcal{K}$. Then we are done if we are able to show $\mathcal{K}=\mathcal{L}$. By the minimality (with respect to set inclusion) of $\mathcal{L}$ among all Lie subalgebras of $\F\monoid/\mathcal{I}_\sigma$ that contain $A,B$, a sufficient condition for $\mathcal{K}$ to be equal to $\mathcal{L}$ is:
\begin{enumerate}[label=\upshape(\roman*), leftmargin=*, widest=iii]
    \item $A,B\in\mathcal{K}$; and\label{hasAB}
    \item $\mathcal{K}$ is a Lie subalgebra of $\F\monoid/\mathcal{I}_\sigma$; and\label{Liesub}
    \item $\mathcal{K}\subseteq\mathcal{L}$.\label{intheLieAlg}
\end{enumerate}
The condition \ref{hasAB} follows immediately from the definition of $\mathcal{K}$. To prove \ref{Liesub}, we show that for any basis vectors $L,R$ of $\mathcal{K}$ among \eqref{LieBasis}, $\left[ L,R\right]$ is a linear combination of \eqref{LieBasis}. We summarize all the cases in the following table.
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|}
\hline
 & $A$ & $B$ & $C^x$ & $C^xA^y$ & $B^yC^x$ \\
 \hline
$A$ & $0$ & $C$ & \eqref{table1} & \eqref{table2} & \eqref{tablei} \\
 \hline
$B$ &  & $0$ & \eqref{table5} & \eqref{tableii} & \eqref{table6} \\
 \hline
$C^k$ &  &  & $0$ & \eqref{table3} & \eqref{table7} \\
 \hline
$C^kA^l$ &  &  &  & \eqref{table4} & Lemma~\ref{tableLem}\\
 \hline
$B^lC^k$ &  &  &  & & \eqref{table8}\\
 \hline
\end{tabular}
\end{center}
All the possiblities for $L$ are written in the first column of the above table, while all the possibilities for $R$ are listed in the first row. The cells that are neither in the first column nor in the first row indicate something about $\left[ L,R\right]$. If the cell indicates $0$ or $C$, then $\left[ L,R\right]$ is equal to the content of that cell. If the cell contains an equation number, then that equation shows how $\left[ L,R\right]$ is a linear combination of \eqref{LieBasis}. The most difficult case corresponds to the cell in the table that indicates ``Lemma~\ref{tableLem},'' and this lemma explains why $\left[ L,R\right]$ for this case is a linear combination of \eqref{LieBasis}. The cells below the main diagonal are left empty because of the skew-symmetry of the Lie bracket. This completes our proof of \ref{Liesub}. To prove \ref{intheLieAlg}, we show that for any basis vector $U$ of $\mathcal{K}$ among \eqref{LieBasis}, $U\in\mathcal{L}$. The proof is clear for the cases $U=A$, $U=B$, $U=C=\left[ A,B\right]$. By Lemma~\ref{CpowerLem}, it can be shown that $C^h\in\mathcal{L}$ for any positive integer $h$ using induction. We are left with the last two possible forms of $U$. By this we mean either $U=C^kA^l$ or $U=B^lC^k$ for some positive integers $k,l$. For these cases, observe that by \eqref{reorderAC}, \eqref{reorderBC}, it can be shown that the relations
\begin{eqnarray}
({\rm ad}\  A)^l(C^k) & = & (m^k-1)^l C^kA^l,\label{adApowerC}\\
({\rm ad}\  B)^l(C^k) & = & (1-m^k)^l B^lC^k, \label{adBpowerC}
\end{eqnarray}
hold in $\F\monoid/\mathcal{I}_\sigma$ for any positive integers $k,l$, using induction. Since we have already established that $C^k\in\mathcal{L}$, the left-hand side of \eqref{adApowerC} is an element of $\mathcal{L}$. Since $m$ is not a root of unity, we can isolate $C^kA^l$ in \eqref{adApowerC}, and conclude that $C^kA^l\in\mathcal{L}$. By similar arguments and computations, the condition $B^lC^k\in\mathcal{L}$ can be deduced from \eqref{adBpowerC}. This proves \ref{intheLieAlg}. Therefore, $\mathcal{K}=\mathcal{L}$.
\end{proof}

\subsection*{Acknowledgements}
This paper was completed while the author was a post-doctoral research guest in M\"alardalen University, V\"aster\aa s, Sweden, hosted by Prof. Sergei Silvestrov, to  whom the author expresses his thanks, and also to Dr. Lars Hellstr\"om for very  valuable insights. The author acknowledges the unerring support of his home department, the Mathematics and Statistics Department of the College of Science, De La Salle University, Manila. This work was partially supported by a grant from the Commission for Developing Countries of the International Mathematical Union (IMU-CDC).


\subsection{Mathematical framework of the forward problem}
\label{subsec:tMatRecovery}
If we define the electric field of light incident on the distal facet of the fibre bundle at a point $(u,v)$ as $T_i (u,v)$ and the electric field incident on the camera sensor at a point $(x,y)$ as $T_o (x,y)$, then we can model the relationship between the two using a spatially dependent point-spread function (or \emph{Green's function}), $h(u,v,x,y)$ as:

\begin{equation}
\label{eq:greenFunc}
T_o(x,y) = \iint h(u, v, x, y) T_i(u,v) \mathrm{d}u \mathrm{d}v
\end{equation}

In the case of free-space propagation, $h(u,v,x,y)$ can be set as $\frac{e^{ikr}}{i r\lambda}$ to give the Fresnel diffraction integral.  However, in our case $h(u,v,x,y)$ is a more complex function that is not known \emph{a priori}. A simple approach to approximate this is to discretise the inputs and outputs and model $h(u,v,x,y)$ with a matrix. To achieve this, we let $T_i(u,v)$ be discretised and then lexicographically ordered into an $n \times 1$ column vector, $\mathbf{x}_h \in \mathbb{C}^n$. We then let $T_o(x,y)$ be discretised and lexicographically ordered into an $m \times 1$ column vector, $\mathbf{y}_h \in \mathbb{C}^m$. In this latter case, an obvious choice of sampling scheme would be to use the values recorded by each pixel of the image sensor. Discretising $h(u,v,x,y)$ leads to the matrix equation:

\begin{equation}
\mathbf{y}_h = \mathbf{A}_{hh} \mathbf{x}_h
\end{equation}

\noindent where $\mathbf{A}_{hh}$ is an $m \times n$ transfer matrix mapping horizontally polarised light at the input to horizontally polarised light at the output.  Repeating the discretisation for the vertically polarised input and output fields, then combining we find:

\begin{equation}
\mathbf{y} =  \left[ \begin{array}{c}
								\mathbf{y}_h \\
								\mathbf{y}_v \end{array} \right]
\end{equation}

\begin{equation}
\mathbf{x} =  \left[ \begin{array}{c}
								\mathbf{x}_h \\
								\mathbf{x}_v \end{array} \right]
\end{equation}

\begin{equation}
\mathbf{A} =  \left[ \begin{array}{cc}
								\mathbf{A}_{hh} & \mathbf{A}_{hv}\\
								\mathbf{A}_{vh} & \mathbf{A}_{vv} \end{array} \right]
\end{equation}

\begin{equation}
\label{eq:forwardEq}
\mathbf{y} = \mathbf{A} \mathbf{x}
\end{equation}

\noindent where $\mathbf{A}$ is the full $2m \times 2n$ transfer matrix.  In practice we record images of the proximal facet of the MCF in the plane of the camera ($\mathbf{y}$) and use these to reconstruct the input fields in the distal facet plane ($\mathbf{x}$).  Equation \ref{eq:forwardEq} can therefore be considered a linear inverse problem: for a known $\mathbf{y}$ and $\mathbf{A}$, the goal is to find $\mathbf{x}$.  One method of solving this problem is to find a pseudo-inverse of $\mathbf{A}$ such that:

\begin{equation}
\mathbf{x} = \mathbf{A}^{-1} \mathbf{y}
\end{equation}

However, we do not typically know $\mathbf{A}^{-1}$ in advance and so we must first determine this before $\mathbf{y}$ can be recovered.  Consider now combining $p$ vectors representing different input fields at the distal facet to form a matrix, $\mathbf{X}$:

\begin{equation}
\label{eq:xMat}
\mathbf{X} =  \left[ \begin{array}{cccc}
\mathbf{x}_1 & \mathbf{x}_2 & \cdots &\mathbf{x}_p
 \end{array} \right]
\end{equation}

For each of these $p$ vectors, the field exiting the distal facet is recorded at the camera plane, creating another matrix, $\mathbf{Y}$:

\begin{equation}
\label{eq:yMat}
\mathbf{Y} =  \left[ \begin{array}{cccc}
\mathbf{y}_1 & \mathbf{y}_2 & \cdots &\mathbf{y}_p
 \end{array} \right]
\end{equation}

Collectively, $\mathbf{X}$ and $\mathbf{Y}$ are termed the \emph{calibration measurements}.  We can then write the following:

\begin{equation}
\mathbf{X} = \mathbf{A}^{-1} \mathbf{Y}
\end{equation}

\begin{equation}
\label{eq:recoveryEq}
\mathbf{X}^T = \mathbf{Y}^T (\mathbf{A}^{-1})^T
\end{equation}

By solving for each of the $p$ columns of $(\mathbf{A}^{-1})^T$ using the respective column of $\mathbf{X}^T$, we can reduce the problem of finding the inverse transmission matrix to $p$ linear inverse problems. Solving each of these provides a single row of $\mathbf{A}^{-1}$:

\begin{equation}
\label{eq:recoveryEq_sing}
\mathbf{x}_r^T = \mathbf{Y}^T \mathbf{a}_r^T
\end{equation}

\noindent where $\mathbf{x}_r$ is the $r^{th}$ row of $\mathbf{X}$ and $\mathbf{a}_r$ is the $r^{th}$ row of $\mathbf{A}^{-1}$.

\subsection{Experimental measurements}
\label{subsubsec:parallelisation}
The full process of producing appropriate sets of inputs and outputs, $\mathbf{X}$ and $\mathbf{Y}$, for recovering $\mathbf{A}^{-1}$ is detailed in Supplementary Figure \ref{fig:matrixRecoveryAlg}. We control the vectors $\mathbf{x}_i$ of Equation \ref{eq:xMat} for $i=1\cdots p$ by displaying an appropriate hologram on SLM1. In this case, each $\mathbf{x}_i$ represents one of 81 translations (over a $9\times9$ grid) in 9$\mu$m steps (limiting spatial resolution) of an array of 12 equispaced spots (Supplementary Figure \ref{fig:holos}a) in one of 3 elliptical polarisation states, resulting in a total of $p=81 \times 3 = 243$ inputs.  Elliptical states are used so as to maximise power transmitted through the fibre, thereby minimising the effect of noise compared to linear states with only half the power.  A minimum of two states are required to resolve the 2 orthogonal polarisation components but further states can reduce noise.  Using 3 states produced a $\sim 33\%$ reduction in noise, but increasing to 4 states provided little improvement as the noise performance was observed to be limited by approximation noise in the basis pursuit algorithm and non-flatness of the hologram.  3 states are therefore selected as a balance between noise reduction and experimental time reduction.  The translation step size is chosen so that every second fibrelet in each dimension is utilised.  This offers a balance between reducing experimental time ($<$1 hour) while still offering sufficient resolution for examining small lesions (size $>50\mu$m).

The choice of 12 spots ensures individual spots are spaced sufficiently far from each other that there is negligible interference between them at the proximal facet of the fibre bundle (even after been scrambled by the matrix $\mathbf{A}$) and they can be considered to behave independently. In a pixel basis representation, the resulting fibre bundle transmission matrix is both \emph{sparse}, and in some sense \emph{local} (technically, a \emph{band matrix}). Hence, not only are most entries zero, but the non-zero entries are clustered together. Consequently, a single image containing 12 spots can be split up to produce 12 separate characterisation images. Each image simply uses the relevant subsection of the full image and sets all other pixels to zero. This means that there are effectively $12 \times 243 = 2916$ characterisation images from only 243 actual measurements. Each of the 81 spatial translations effectively provides 12 pixels of resolution at the distal facet (one per spot).  Therefore, in total there are $12 \times 81 = 972$ pixels that can be spatially resolved at the distal facet.

When using a phase retrieval algorthm (see Online Methods), the phase recovered for each image has an arbitrary global offset. If measured na\"{\i}vely, the columns of $\mathbf{Y}$ in Equation \ref{eq:yMat} would have an arbitrary phase offset between them. To compare different phase images, a fixed phase reference is created during the characterisation measurements by displaying a superposition of two holograms on SLM1: the first hologram creates the translatable array of spots and the second creates a single spot that is kept in a fixed position. The fixed position is a set of fibrelets at the edge of the MCF, as shown in Supplementary Figure \ref{fig:matrixRecoveryAlg}a. Recovered phase values are normalised to the phase of this spot at the proximal facet of the fibre. This is effectively a form of common-path interferometry. Such a global phase offset is also termed \emph{piston} and can be thought of as a phase surface polynomial of order 0 \cite{Loterie2015a}. Here, we observed experimentally that in addition to this 0$^{th}$ order correction, a 1$^{st}$ order time-dependent phase correction,termed \emph{phase tilt}, was also required. The origin of the time-dependent phase tilt is likely due to temperature, affecting the SLM and optical mounts \cite{Garcia-Marquez2011}. Air currents could also play a role, however the experiment is entirely contained in an optical box so their effect should be negligible. The 1$^{st}$ order phase correction term is tracked by measuring a reference image at 60 second intervals during characterisation. The relative phase tilt between reference images is computed and interpolation with respect to time is then used to estimate the correction to be applied to characterisation measurements, $\mathbf{y_i}$.

\subsection{Solving the inverse problem}
The region of the image sensor used for data collection is 1200$\times$1200 pixels and so the inverse transmission matrix, $\mathbf{A}^{-1}$, has 2$\times$1200$\times$1200 =  2.88 million columns when the two polarisation states are considered. Similarly, the 81 translations of 12 spots in three polarisation states at the proximal facet gives a total of 2916 rows for $\mathbf{A}^{-1}$. To store this 2916 $\times$2880000 matrix at single floating point precision would require $>$30Gb of memory. This can make recovery of the matrix and subsequent operations infeasibly slow.

However, because the transmission matrix is sparse and the few non-zero elements are grouped together (i.e. a banded structure as discussed above), it is possible to reconstruct only the parts of $\mathbf{A}^{-1}$ where the non-zero elements are located. Consider again Equation \ref{eq:recoveryEq_sing} in which we seek to solve for row $r$ of $\mathbf{A}^{-1}$, $\mathbf{a}_r$:

\begin{equation}
\mathbf{x}_r^T = \mathbf{Y}^T \mathbf{a}_r^T
\end{equation}

Here, $\mathbf{x}_r$ is a single row of $\mathbf{X}$ and thus corresponds to a single spatial position across all input calibration fields at the distal facet.  Similarly, $\mathbf{a}_r$ is a single row of $\mathbf{A}^{-1}$ and therefore represents the coupling to the single spatial position represented by $\mathbf{x}_r$ from all spatial positions at the proximal facet, represented columns of $\mathbf{Y}$.  However, because we know $\mathbf{A}^{-1}$ is sparse and banded, the only contributions to $\mathbf{x}_r$ will come from points in a confined region on the distal facet.  The distal and proximal facet coordinate axes are the same as a result of an affine transform (see Online Methods) and so if the coordinates of the spatial position of $\mathbf{x}_r$ are $(a,b)$, then the relevant points on the distal facet will be neighbouring $(a,b)$.  We therefore select a set of spatial points on the distal facet (i.e. columns $\mathbf{Y}$) that are located within a square centred on $(a,b)$ with side length 68.2$\mu$m.  This length is exactly equal to the spacing between adjacent spots in the array used to characterise the fibre (Supplementary Figure \ref{fig:holos}a), which was in turn set so as to avoid interference of spots at the distal facet (Online Methods). This process is illustrated in Supplementary Figure \ref{fig:matrixRecoveryAlg}.  Equation \ref{eq:recoveryEq_sing} then becomes:

\begin{equation}
\label{eq:subEq}
\mathbf{x}_{r}^T = \mathbf{Y}_{sub}^T \mathbf{a}_{r,sub}^T
\end{equation}

\noindent where $\mathbf{Y}_{sub}$ is the selected subset of rows of $\mathbf{Y}$ and
 $\mathbf{a}_{r,sub}$ is the corresponding subset of columns of $\mathbf{A}^{-1}$. Typically, $\mathbf{a}_{r,sub}$ has around 60,000 columns over 2 polarisations, a 50-fold reduction from the original 2.88 million.  However, given that $\mathbf{x}_r$ has only $\sim$3000 columns, the problem is underdetermined by a factor of 10.  The problem can still be solved by exploiting the prior knowledge that each $\mathbf{a}_{r,sub}$ ought to be sparse because spatial positions towards the edge would likely couple little or no power directly from the central positions (in agreement with observations). This is achieved using an $L1$ regulariser so that the solution to Equation \ref{eq:subEq} becomes \cite{Zhang2015}:

\begin{equation}
\min ||\mathbf{a}_{r,sub}||_1 \mathrm{~~subject~to~~} ||\mathbf{x}_r^T - \mathbf{Y}_{sub}^T \mathbf{a}_{r,sub}^T||_2 < \delta
\label{eq:solver}
\end{equation}

\noindent where $\delta$ is a parameter that can be adjusted to favour a more or less sparse solution.  Equation \ref{eq:solver} is solved numerically using a basis pursuit denoising algorithm implemented in the SPGL1 package \cite{spgl1:2007}.

\pagebreak
\noindent
\Large\textbf{Supplementary Note 2:}\\
\Large\textbf{Evaluation of polarisation parameters}\\
\normalsize
\setcounter{section}{2}
\setcounter{subsection}{0}
\subsection{Polarisation model}
\label{subsec:polModel}
When polarised light interacts with the turbid medium of tissue, it may simultaneously encounter: depolarisation, arising from multiple scattering of sub-cellular structures; birefringence and diattenuation, arising from anisotropic muscle fibres and structural proteins; and circularity caused by chiral molecules and structures.

The propagation of polarised light and its interaction within a sample can be described by either Jones (field-based) or Stokes-Mueller (intensity-based) calculus \cite{Simon1982}. Stokes-Mueller is more commonly used for tissue polarimetry as it can handle partial polarisation states, represented by the depolarisation parameter  \cite{Ghosh2011,Qi2016,Yao1999,Pierangelo2013a}, though Jones calculus has found use in polarisation-sensitive optical coherence tomography of tissue \cite{Park2009,Li2015}.  However, if the coherence length of the light source is substantially longer than possible scattering paths in the tissue, then the light always has a well-defined polarisation state \cite{Simon1982}.  This corresponds to a light source with a high degree of temporal coherence, i.e. monochromatic light. The laser used here has a coherence length of around 10 m compared to the typical transport mean free path for tissue of $<$ 1 mm \cite{Ntziachristos2010} and so this condition is satisfied. In addition to this temporal coherence, the light field must also exhibit spatial coherence in order to have a well-defined polarisation state.  This can be violated, for example, if the integration area of each camera pixel contains  different polarisation states.  However, in our case the light field from the sample is first passed through the fibre bundle where the single-mode propagation in each fibrelet admits only a single polarisation state -- a filter of sorts.  Each fibrelet is imaged on the camera using $>$10 pixels so it can therefore be assumed that the polarisation state is uniform within each pixel area. Depolarisation is therefore considered to be negligible in this work.

In order to examine the properties of birefringence, diattenuation (also called \emph{linear dichroism}) and circularity, we choose a generalised decomposition of the Jones matrix into an arbitrary elliptical retarder followed by a partial linear polariser. The elliptical retarder accounts for birefringent retardance and circularity. The partial linear polariser accounts for the fact that anisotropic molecules might absorb light preferentially along some axis. This is illustrated in Supplementary Figure \ref{fig:polModel}. It is noted that this model does not account for some polarimetric properties, notably circular dichroism. Circular dichroism has applications in spectroscopic analysis of chiral molecules (e.g. many proteins and biological molecules) but produces weak signals so is typically only effective for molecules in concentrated solution \cite{Fasman1996}. Further, little or no work has been produced to date linking circular dichroism at optical wavelengths to imaging of cancer. Therefore, for simplicity it is neglected in this work.

Having measured a set of $n$ output polarisation vectors from the sample $\mathbf{V}(x,y) = \left[\mathbf{v}_1(x,y) \cdots \mathbf{v}_n(x,y)\right]$ and knowing the corresponding $n$ sample illumination polarisation vectors $\mathbf{U}(x,y) = \left[\mathbf{u}_1(x,y) \cdots \mathbf{u}_n(x,y)\right]$ the 2$\times$2 complex Jones matrices $\mathbf{J}(x,y)$ representing how the sample maps input polarisation state to output polarisation state at each pixel can be determined through:

\begin{equation}
\label{eq:forwardModelJones}
\mathbf{V}(x,y) = \mathbf{J}(x,y) \mathbf{U}(x,y)
\end{equation}

\noindent where $\mathbf{J}(x,y)$ is the Jones matrix of a particular pixel at $(x,y)$. We express our decomposition of $\mathbf{J}(x,y)$ as:

\begin{equation}
\label{eq:jonesMat}
\mathbf{J}(x,y) = \mathbf{A}_{pol} \mathbf{A}_{ret}
\end{equation}

\begin{equation}
\label{eq:polMat}
\mathbf{A}_{pol} = \left( \begin{array}{cc}
\cos \theta_D & \sin \theta_D \\
-\sin \theta_D & \cos \theta_D \end{array} \right)
\left( \begin{array}{cc}
\sqrt{1+D} & 0 \\
0 & \sqrt{1-D} \end{array} \right)
\left( \begin{array}{cc}
\cos \theta_D & -\sin \theta_D \\
\sin \theta_D & \cos \theta_D \end{array} \right)
\end{equation}

\begin{equation}
\label{eq:retMat}
\mathbf{A}_{ret} = \left( \begin{array}{cc}
1 & 0 \\
0 & e^{i\xi}\end{array} \right)
\left( \begin{array}{cc}
\cos \theta_{\eta} & \sin \theta_{\eta} \\
-\sin \theta_{\eta} & \cos \theta_{\eta} \end{array} \right)
\left( \begin{array}{cc}
e^{i\eta/2} & 0 \\
0 & e^{-i\eta/2} \end{array} \right) \cdot \left( \begin{array}{cc}
\cos \theta_{\eta} & -\sin \theta_{\eta} \\ \sin \theta_{\eta} & \cos \theta_{\eta} \end{array} \right)
\left( \begin{array}{cc}
1 & 0 \\
0 & e^{-i\xi}\end{array} \right)
\end{equation}

\noindent allowing us to recover the following 5 polarisation properties: diattenuation, $D$, diattenuation axis orientation, $\theta_D$, retarder circularity, $\xi$, retardance, $\eta$, and retardance axis orientation, $\theta_{\eta}$.  From these equations, it is observed that the inferred parameters are all real numbers limited to the following ranges:

\begin{equation}
\label{eq:paramRanges}
\theta_{\eta} \in (-\pi/2, \pi/2],~~~~~~~\eta \in (-\pi, \pi],~~~~~~~\xi \in (-\pi, \pi], ~~~~~~~D \in [-1, 1],~~~~~~~\theta_D \in (-\pi/2, \pi/2]
\end{equation}

\subsection{Ordering of polarisation elements}
\label{subsec:polOrdering}
The model for inference of polarimetric parameters, introduced in Section \ref{subsec:polModel} and illustrated in Supplementary Figure \ref{fig:polModel}, comprises an elliptical retarded followed by a partial polariser. In the Jones formalism these are each represented by $2\times2$ matrices, which are multiplied together to obtain overall behaviour.  However, matrix multiplication is non-commutative and so it is necessary to consider an alternative model in which the order of these two elements is reversed.  By adapting Equation \ref{eq:jonesMat} we get an alternative composite Jones matrix:

\begin{equation}
\label{eq:jonesMat_rev}
\mathbf{\hat{J}}(x,y) = \mathbf{A}_{ret} \mathbf{A}_{pol}
\end{equation}

If we then expand both Equation \ref{eq:jonesMat} and \ref{eq:jonesMat_rev} using Equations \ref{eq:polMat} and \ref{eq:retMat}, we get:

\begin{equation}
\label{eq:polOrder1}
\mathbf{J}(x,y) = \left(\begin{array}{c c}
K & L \\
M & N
\end{array}\right)
\end{equation}

\noindent and

\begin{equation}
\mathbf{\hat{J}}(x,y) = \left(\begin{array}{c c}
\hat{K} & \hat{L} \\
\hat{M} & \hat{N}
\end{array}\right)
\end{equation}

\noindent where

\begin{equation}
\begin{array}{c}
K = \left({\sin\left(\theta_D\right)}^2 \sqrt{1 - D} + {\cos\left(\theta_D\right)}^2 \sqrt{D + 1}\right) \left(\cos\left(\frac{\eta}{2}\right) + \cos\left(2 \theta_{\eta}\right) \sin\left(\frac{\eta}{2}\right) i\right) \\
- \frac{1}{2} \sin\left(\frac{\eta}{2}\right) \sin\left(2 \theta_{\eta}\right) \sin\left(2\theta_D\right) \left(\sqrt{D + 1} - \sqrt{1 - D}\right) \left(\sin\left(\xi\right) + i\cos\left(\xi\right)\right)
\end{array}
\end{equation}

\begin{equation}
\begin{array}{c}
\hat{K} = \left({\sin\left(\theta_D\right)}^2 \sqrt{1 - D} + {\cos\left(\theta_D\right)}^2 \sqrt{D + 1}\right) \left(\cos\left(\frac{\eta}{2}\right) + \cos\left(2 \theta_{\eta}\right) \sin\left(\frac{\eta}{2}\right) i\right) \\
- \frac{1}{2} \sin\left(\frac{\eta}{2}\right) \sin\left(2 \theta_{\eta}\right) \sin\left(2\theta_D\right) \left(\sqrt{D + 1} - \sqrt{1 - D}\right) \left(- \sin\left(\xi\right) + i\cos\left(\xi\right)\right)
\end{array}
\end{equation}

\begin{equation}
\begin{array}{c}
L = i \mathrm{e}^{i \xi} \sin\left(\frac{\eta}{2}\right) \sin\left(2\theta_{\eta}\right) \left({\sin\left(\theta_D\right)}^2 \sqrt{1 - D} + {\cos\left(\theta_D\right)}^2 \sqrt{D + 1}\right) \\
- \frac{1}{2} \sin\left(2\theta_D\right) \left(\sqrt{D + 1} - \sqrt{1 - D}\right) \left(\cos\left(\frac{\eta}{2}\right) - i\cos\left(2 \theta_{\eta}\right) \sin\left(\frac{\eta}{2}\right)\right)
\end{array}
\end{equation}

\begin{equation}
\begin{array}{c}
\hat{L} = i\mathrm{e}^{i\xi} \sin\left(\frac{\eta}{2}\right) \sin\left(2 \theta_{\eta}\right) \left({\sin\left(\theta_D\right)}^2 \sqrt{D + 1} + {\cos\left(\theta_D\right)}^2 \sqrt{1 - D}\right) \\
- \frac{1}{2}\sin\left(2 \theta_D\right) \left(\sqrt{D + 1} - \sqrt{1 - D}\right) \left(\cos\left(\frac{\eta}{2}\right) + i\cos\left(2 \theta_{\eta}\right) \sin\left(\frac{\eta}{2}\right)\right)
\end{array}
\end{equation}

\begin{equation}
\begin{array}{c}
M = \mathrm{e}^{i \xi} \sin\left(\frac{\eta}{2}\right) \sin\left(2 \theta_{\eta}\right) \left({\sin\left(\theta_D\right)}^2 \sqrt{D + 1} + {\cos\left(\theta_D\right)}^2 \sqrt{1 - D}\right) \\
- \cos\left(\theta_D\right) \sin\left(\theta_D\right) \left(\sqrt{D + 1} - \sqrt{1 - D}\right) \left(\cos\left(\frac{\eta}{2}\right) + i\cos\left(2 \theta_{\eta}\right) \sin\left(\frac{\eta}{2}\right)\right)
\end{array}
\end{equation}

\begin{equation}
\begin{array}{c}
\hat{M} = \mathrm{e}^{i \xi} \sin\left(\frac{\eta}{2}\right) \sin\left(2 \theta_{\eta}\right) \left({\sin\left(\theta_D\right)}^2 \sqrt{1 - D} + {\cos\left(\theta_D\right)}^2 \sqrt{D + 1}\right) \\
-  \cos\left(\theta_D\right)\sin\left(\theta_D\right) \left(\sqrt{D + 1} - \sqrt{1 - D}\right) \left(\cos\left(\frac{\eta}{2}\right) + i\cos\left(2 \theta_{\eta}\right) \sin\left(\frac{\eta}{2}\right)\right)
\end{array}
\end{equation}

\begin{equation}
\begin{array}{c}
N = \left({\sin\left(\theta_D\right)}^2 \sqrt{D + 1} + {\cos\left(\theta_D\right)}^2 \sqrt{1 - D}\right) \left(\cos\left(\frac{\eta}{2}\right) + i\cos\left(2 \theta_{\eta}\right) \sin\left(\frac{\eta}{2}\right)\right)\\
 -\frac{1}{2}\sin(\frac{\eta}{2})\sin\left(2\theta_D\right)  \sin\left(2\theta_{\eta}\right)\left( \sqrt{D + 1} - \sqrt{1 - D}\right) \left(-\sin(\xi)+i\cos(\xi) \right)
\end{array}
\end{equation}

\begin{equation}
\begin{array}{c}
\hat{N} = \left({\sin\left(\theta_D\right)}^2 \sqrt{D + 1} + {\cos\left(\theta_D\right)}^2 \sqrt{1 - D}\right) \left(\cos\left(\frac{\eta}{2}\right) + i\cos\left(2 \theta_{\eta}\right) \sin\left(\frac{\eta}{2}\right)\right)\\
 - \frac{1}{2}\sin\left(\frac{\eta}{2}\right) \sin\left(2 \theta_D\right) \sin\left(2 \theta_{\eta}\right) \left(\sqrt{D + 1} - \sqrt{1 - D}\right) \left(\sin\left(\xi\right) + i\cos\left(\xi\right)\right)
\end{array}
\end{equation}

Taking the difference between these two matrices, we find that:

\begin{equation}
\label{eq:polDiffMat}
\mathbf{J}(x,y) - \mathbf{\hat{J}}(x,y) = \mathbf{J'}(x,y) = \left(\begin{array}{c c}
K' & L' \\
M' & N'
\end{array}\right)
\end{equation}

\noindent where

\begin{equation}
\begin{array}{c}
K' = K - \hat{K} = - \left(\sqrt{D + 1} - \sqrt{1 - D}\right) \sin\left(\frac{\eta}{2}\right) \sin\left(2 \theta_{\eta}\right) \sin\left(2\theta_D\right) \sin\left(\xi\right)
\end{array}
\end{equation}

\begin{equation}
\begin{array}{c}
L' = L - \hat{L} = \left(\sqrt{D + 1} - \sqrt{1 - D}\right)\left( i \mathrm{e}^{i \xi} \sin\left(\frac{\eta}{2}\right) \sin\left(2\theta_{\eta}\right)\cos(2\theta_D)+ i\sin\left(2\theta_D\right) \left(\cos\left(2 \theta_{\eta}\right) \sin\left(\frac{\eta}{2}\right)\right)\right)
\end{array}
\end{equation}

\begin{equation}
\begin{array}{c}
M' = M - \hat{M} = -\left(\sqrt{D + 1} - \sqrt{1 - D}\right) \mathrm{e}^{i \xi} \sin\left(\frac{\eta}{2}\right) \sin\left(2 \theta_{\eta}\right) \cos\left(2\theta_D\right)
\end{array}
\end{equation}

\begin{equation}
\begin{array}{c}
N' = N - \hat{N} = \left( \sqrt{D + 1} - \sqrt{1 - D}\right) \sin(\frac{\eta}{2})\sin\left(2\theta_D\right)  \sin\left(2\theta_{\eta}\right) \sin(\xi)
\end{array}
\end{equation}

The error terms $K'$, $L'$, $M'$ and $N'$ are all proportional to $\left( \sqrt{D + 1} - \sqrt{1 - D}\right)$, and so as $D \to 0$, so these error terms become negligible.  To investigate the range of $D$ for which the error between these two models is negligible, we performed a Monte-Carlo simulation over the broadest possible range of parameters:

\begin{equation*}
\theta_{\eta} \sim U(-\pi/2,\pi/2)
\end{equation*}

\begin{equation*}
\eta \sim U(-\pi,\pi)
\end{equation*}

\begin{equation*}
\xi \sim U(-\pi,\pi)
\end{equation*}

\begin{equation*}
\theta_D \sim U(-\pi/2,\pi/2)
\end{equation*}

\begin{equation*}
D \sim N(0,\sigma_D)
\end{equation*}

100,000 samples were drawn from these distributions for different values of $\sigma_D$ using a MATLAB script.  For each sample drawn, Jones matrices are computed for the two orderings ($\mathbf{J}$ and $\mathbf{\hat{J}}$ respectively) of the polariser and retarder, enabling computation of the error matrix $\mathbf{J}'$ of Equation \ref{eq:polDiffMat}.  The elements are then normalised and summed giving an error metric, $Q$, defined (with reference to Equations \ref{eq:polOrder1} and \ref{eq:polDiffMat}) as:

\begin{equation*}
Q = \frac{1}{4} \left( \left|\frac{K'}{K}\right| + \left|\frac{L'}{L}\right| + \left|\frac{M'}{M}\right| + \left|\frac{N'}{N}\right| \right)
\end{equation*}

\begin{equation}
Q = \frac{1}{4} \left( \left|\frac{K - \hat{K}}{K}\right| + \left|\frac{L - \hat{L}}{L}\right| + \left|\frac{M - \hat{M}}{M}\right| + \left|\frac{N - \hat{N}}{N}\right| \right)
\end{equation}

This represents the average absolute proportional deviation introduced by reversing polarisation elements.  The measured cumulative distribution function for different values of $\sigma_D$ is shown in Supplementary Figure \ref{fig:polMonteCarlo}.  In the samples of mouse oesophageal tissue presented here, the true value of $\sigma_D$ is estimated to be 0.0081.  As this is less than 0.01, it is observed from the cumulative density function that in $>$99.5\% of possible cases the error introduced by a different ordering of polarisation elements will be $<$5\%. If a 10\% error threshold is considered acceptable, then a value of $\sigma_D = $0.05, 5-times the measured value, would still be sufficiently accurate in $>$96\% of cases.  Of course, if samples are expected to have larger values of $D$ then the two orderings would give different results and it would be necessary to determine which ordering had the stronger biological underpinning and best explained the measured data.  Such a `model comparison' could be achieved through an extension of the Bayesian inference approach of Supplementary Section \ref{subsec:bayesInf}.

\subsection{Bayesian inference}
\label{subsec:bayesInf}
The polarisation model is fitted to the measured data in order to extract the 5 polarisation properties ($\theta_D$, $D$, $\xi$, $\theta_{\eta}$, $\eta$) using a Bayesian inference approach, which has been found to be more robust to noise than clustering or pseudo-inverse approaches \cite{Zallat2008,Kasaragod2014}. In this work maximum likelihood estimates are determined during this inference and are used to plot images.  These could be computed more quickly by deriving analytical expressions via eigenvalue decomposition of the Jones matrix \cite{Makita2010, Lu1996}.  However, such an approach does not easily permit estimation of joint posterior probabilities of neighbouring pixels, which is essential for compensating beam misalignment (see Section \ref{subsubsec:beamMisalignment}). Therefore, we numerically infer joint maximum likelihood values.

To fit the model, we evaluate the \emph{posterior} probability of the parameter set as a function of the data using Bayes' theorem:

\begin{align*}
p\left[ D, \theta_D, \eta, \theta_{\eta}, \xi | \mathbf{U}(x,y), \mathbf{V}(x,y)\right. & \left.\right] \propto \\
p\left[\mathbf{V}(x,y) = \right. & \left. \mathbf{J}(D, \theta_D, \eta, \theta_{\eta}, \xi,x,y) \mathbf{U}(x,y)\right] \cdot p(D, \theta_D, \eta, \theta_{\eta}, \xi)
\end{align*}

The values in $\mathbf{V}(x,y)$ represent measured complex quantities and so are assumed to be independent and drawn from complex Gaussian distributions as:

\begin{equation}
\label{eq:complexNormal}
\mathbf{v}(x,y) \sim \mathcal{CN}\left[\mathbf{J}(D, \theta_D, \eta, \theta_{\eta}, \xi,x,y) \mathbf{u}(x,y), \sigma^2 \mathbf{I}\right]
\end{equation}

\noindent where $\mathbf{v}(x,y)$ is a column of $\mathbf{V}(x,y)$, $\mathbf{u}(x,y)$ is a column of $\mathbf{U}(x,y)$, $\sigma^2 \mathbf{I}$ is the covariance matrix, and $\mathcal{CN} (\mathbf{\mu},\mathbf{\Sigma})$ is a 2-D complex Gaussian distribution of mean $\mathbf{\mu}$ and covariance $\mathbf{\Sigma}$. The noise standard deviation, $\sigma$ is an additional parameter that can be inferred from the data.  However, after initial processing of several samples, the value was fixed at its inferred maximum likelihood value of 0.4 for subsequent calculations in order to increase speed.  We wish to use a broad prior to minimise any bias in the results.  The broadest possible prior would be a 5-dimensional uniform distribution over all the polarisation parameters (ignoring $\sigma$ by assuming it is now fixed).  In this case, the prior distributions for each individual parameter are independent so that we can write:

\begin{equation}
p(D, \theta_D, \eta, \theta_{\eta}, \xi) = p(D)p(\theta_D)p(\eta)p(\theta_{\eta})p(\xi)
\end{equation}

It is assumed that in the case of a slightly more restrictive prior, this independence assumption will still be valid to a good approximation and so it used for the rest of the analysis presented here.  Priors for parameters representing angles, i.e. circular quantities, are modelled with \emph{von Mises} distributions. The von Mises distribution is a continuous probability distribution defined on a circular domain of $[-\pi,\pi]$ that wraps around, giving a probability density function of $\frac{e^{\kappa \cos (x-\mu)}}{2 \pi I_0(\kappa)}$ where $\mu$ is the mean, $1/\kappa$ is analogous to variance, and $I_0(\cdots)$ is a modified Bessel function of order 0. For $\theta_D$, $\theta_{\eta}$ and $\eta$, the $\kappa$ is set to $\infty$ providing the broadest possible prior. For $\xi$, a restricted value of $\kappa = 1$ is used.  This is because when performing inference of known test targets, it was found that this parameter is prone to over-fitting and that a slightly restrictive prior substantially prevents this.  Further, previous studies that have shown that even in samples with high concentrations of chiral molecules such as glucose, the measured circularity is still relatively small \cite{Guo2007}.  It is therefore not expected that this parameter will be very large for the phantoms or tissue samples measured. Diattenuation, $D$, is the only non-circular parameter (i.e. it is not an angle) and so its prior is a truncated Gaussian over its domain $[-1,1]$ with mean 0 and broad variance.  Again, the broad (as opposed to infinite) variance is used to prevent overfitting and is consistent with previous work that has found the values of $D$ to be small (typically $<$0.2) in real tissue samples \cite{Qi2016,Pierangelo2013a}.  It is noted that the ability to find more accurate solutions in the presence of noise or incomplete information using such prior knowledge is a major benefit of the Bayesian inference approach.

The inference of parameters from data was performed numerically using the STAN package \cite{StanMan}, which produces samples from \emph{a posteriori} distributions using a Markov Chain Monte Carlo simulation engine.  It can also provide maximum likelihood estimates through an optimisation engine, which is used here on account of its high speed, and can also perform variational inference.

\subsection{Compensation for beam misalignment}
\label{subsubsec:beamMisalignment}
The holographic endoscope design exploits two orthogonally polarised arms that are split then recombined before illuminating the sample. This is achieved using a custom housing for the mirrors and polarising beam splitters, which provides a compact way of transferring the beam path between two optical axes at different heights. As a result of this dual arm approach, any slight misalignment of the two polarisation arms leads to a spatially varying illumination polarisation state at the sample, illustrated in Supplementary Figure \ref{fig:waveplateImaging}a. For samples with small retardance this misalignment is not apparent because it is present for both the characterisation of the fibre transmission matrix and the sample imaging steps, and so the relative change in polarisation state is correct.

However, for samples that exhibit large values of retardance, $\eta$, and in which the retardance axes are not aligned with the polarisation axes of the illumination (i.e. ($\theta_{\eta} \neq 0,\pm\pi/2$), light is coupled between the slightly misaligned horizontally and vertically polarised beams. This produces, in the illumination polarisation basis, wavefronts with different tilts in the two polarisation arms.  This manifests as an apparent spatially varying retardance, illustrated in Supplementary Figure \ref{fig:waveplateImaging}a. This apparent retardance does not appear when the polarisation axes are aligned with the retardance axes. Therefore, the effect is removed by re-expressing Equation \ref{eq:forwardModelJones} in a linear polarisation basis at angle $\hat{\theta}_{\eta}$ to the illumination polarisation axes:

\begin{equation}
\label{eq:rotatedPolAxes}
\mathbf{V}(x,y) = \left. \mathbf{J}(x,y) \right|_{\theta_{\eta}=0}
\left( \begin{array}{cc}
\cos \hat{\theta}_{\eta} & \sin \hat{\theta}_{\eta} \\
-\sin \hat{\theta}_{\eta} & \cos \hat{\theta}_{\eta} \end{array} \right)
\mathbf{U}(x,y)
\end{equation}

We then apply Bayesian inference in this new, rotated basis with $\hat{\theta}_{\eta}$ taking the place of $\theta_{\eta}$ as the inferred orientation angle of the retardance axes. To find the $\hat{\theta}_{\eta}$ that minimises the spatial tilt of the retardance (confirming the new basis is aligned with the retardance axes) we compute the joint probability for a small set of $R$ neighbouring pixels surrounding $(x,y)$.  For computational tractability, we employ a technique from variational Bayesian inference that approximates the posterior probability density function as a product of distributions \cite{MacKay2003}:

\begin{align*}
p[ D(x_1,y_1)\cdots D(x_R, y_R), &\theta_D(x_1,y_1)\cdots, \eta(x_1,y_1)\cdots, \hat{\theta}_{\eta}(x_1,y_1)\cdots, \xi(x_1,y_1)\cdots | \mathbf{U}(x_1,y_1)\cdots, \mathbf{V}(x_1,y_1)\cdots] \\
&= \prod_{r=1}^R p[D(x_r,y_r),\theta_D(x_r,y_r),\eta(x_r,y_r), \hat{\theta}_{\eta}(x_r,y_r), \xi(x_r,y_r) | \mathbf{U}(x_r,y_r), \mathbf{V}(x_r,y_r)]
\end{align*}

\noindent where the term inside the product is computed using equations \ref{eq:complexNormal} and \ref{eq:rotatedPolAxes}. This approximation implies independence of the conditional probabilities of polarisation parameters of neighbouring pixels.  This is a reasonable assumption because the conditional probability distributions (that is, the probability distributions given the experimentally measured Jones matrices) are effectively accounting for noise, which is assumed to be produced by an independent random process at each pixel.

The maximum likelihood probability is maximised when neighbouring pixels have similar conditional distributions for their polarisation parameters. This implicitly biases towards values of $\hat{\theta}_{\eta}$ where the spatial variation of retardance is minimised, i.e. when the means of the retardance distributions across neighbouring pixels are most similar.

It is not known \emph{a priori} whether the phantoms and tissue samples measured will exhibit significant birefringence.  However, the advantage of the above approach is that it will still provide the best estimate of parameters even if no birefringent tilt artefact is present. This is because the probability distribution of $\hat{\theta}_{\eta}$ tends towards being uniform when $\eta \xrightarrow{}$0, which is the same behaviour expected for $\theta_{\eta}$ (i.e. without the correction).  Therefore, the correction is applied to all measurements taken and will have more effect where it is naturally required.

The number of adjacent pixels used for the correction is set as $R$ = 15.  For optimal resolution, this value should be as small as possible because it effectively reduces resolution by low-pass filtering the image.  On the other hand, larger values can correct for smaller tilt artefacts.  $R$ = 15 was found experimentally to be the smallest value that was sensitive enough to correct for the degree of tilt artefact typically observed for the samples measured (which is in turn dependent on the degree of system misalignment). The outcome of this correction is validated by rotating a half-waveplate in the sample plane (Supplementary Figure \ref{fig:waveplateImaging}).

\subsection{Degeneracy of solutions}
We are decomposing a single matrix into the product of 8 matrices (Equations \ref{eq:jonesMat}, \ref{eq:polMat} and \ref{eq:retMat}), each comprising sinusoidal or complex exponential terms, hence there will invariably be degeneracy of the solutions. For example, consider the matrix $\mathbf{A}_{ret}$.  The 5 matrices of equation \ref{eq:retMat} can be multiplied to give:

\begin{equation}
\mathbf{A}_{ret} = \left( \begin{array}{cc}
e^{i\eta/2} \cos^2 \theta_{\eta} + e^{-i\eta/2} \sin^2 \theta_{\eta} & \left(e^{i\eta/2} - e^{-i\eta/2}\right) e^{-i \xi} \cos \theta_{\eta} \sin \theta_{\eta} \\
\left(e^{i\eta/2} - e^{-i\eta/2}\right) e^{i \xi} \cos \theta_{\eta} \sin \theta_{\eta} & e^{i\eta/2} \sin^2 \theta_{\eta} + e^{-i\eta/2} \cos^2 \theta_{\eta}  \end{array} \right)
\end{equation}

By substituting $\theta_{\eta} \rightarrow \theta_{\eta} + \pi/2$ and $\eta \rightarrow -\eta$, the first row, first column element of $\mathbf{A}_{ret}$ becomes:

\begin{equation*}
e^{-i\eta/2} \cos^2 (\theta_{\eta}+\pi/2) + e^{i\eta/2} \sin^2 (\theta_{\eta}+\pi/2)= e^{i\eta/2} \cos^2 \theta_{\eta} + e^{-i\eta/2} \sin^2 \theta_{\eta}
\end{equation*}

The element of $\mathbf{A}_{ret}$ is therefore unchanged; the other 3 elements are also unchanged under this substitution.  Therefore, the parameter sets  $(\theta_{\eta}, \eta)$ and $(\theta_{\eta} + \pi/2, -\eta)$ are degenerate.

Considering the ranges of all parameters, outlined in Equation \ref{eq:paramRanges}, there is in fact an 8-fold ambiguity for each point in the 5-D parameter space, comprising the following parameter sets:

\begin{equation}
\begin{array}{ccccccc}
(&\theta_{\eta}, &\eta,&\xi,&D,&\theta_D&)\\
(&\theta_{\eta}+\pi/2,&-\eta,&\xi,&D,&\theta_D&)\\
(&-\theta_{\eta},&\eta,&\xi+\pi,&D,&\theta_D&)\\
(&-\theta_{\eta}+\pi/2,&-\eta, &\xi+\pi,&D, &\theta_D&)\\
(&\theta_{\eta},&\eta, &\xi, &-D, &\theta_D + \pi/2&)\\
(&\theta_{\eta}+\pi/2, &-\eta, &\xi, &-D, &\theta_D + \pi/2&)\\
(&-\theta_{\eta}, &\eta, &\xi+\pi, &-D, &\theta_D + \pi/2&)\\
(&-\theta_{\eta}+\pi/2, &-\eta, &\xi+\pi, &-D, &\theta_D + \pi/2&)
\end{array}
\end{equation}

\noindent where parameters are defined to wrap around their ranges, listed in Equation \ref{eq:paramRanges}, via a modulus operation. If we were to randomly select one of the 8 degenerate solutions for each pixel, it could give the misleading impression that polarisation properties exhibit a high degree of variation.  To avoid this scenario, we select a point in 5-D parameter space (here we use $(\theta_{\eta} = \pi/4, \eta = \pi, \xi = 0, D = 1, \theta_D = \pi/4)$) and choose the degenerate solution closest to that point.  This ensures the maximum degree of smoothness in the resulting polarisation images. As a result, observed variations can be confidently classed as structural tissue features rather than artefacts.

\pagebreak
\noindent
\Large\textbf{Supplementary Note 3:}\\
\Large\textbf{Preparation and characterisation of tissue mimicking phantoms}\\
\normalsize
\setcounter{section}{3}
\setcounter{subsection}{0}

\subsection{Preparation of optically scattering and birefringent phantoms}
All chemicals were purchased from Sigma Aldrich unless otherwise stated. Birefringent phantoms were made from acrylamide/bis-acrylamide (30\% concentration) solution as the base material. Ammonium persulfate (APS) and tetramethylethylenediamine (TEMED) were used as catalysts for the co-polymerisation of bis-acrylamide and acrylamide gels. Phantoms were fabricated inside a ventilated fumehood cupboard by mixing 30 mL of Acrylamide/Bis-acrylamide with 600 $\mu$L of APS and 40 $\mu$L of TEMED under thorough vortexing.  The mixture was poured immediately into petri dishes and left for approximately 30 minutes until the polymerisation was complete. Birefringent phantoms having thickness $\sim$2 mm were fabricated by pouring 5 mL of the mixture into the petri dishes. The set phantoms were then extracted from the petri dishes and were mounted on to the mechanical stretcher assembly as shown in Supplementary Figure \ref{fig:stretcher}.

Agar-based scattering phantoms were fabricated using 1.5 w/v\% agar solution as the base material as described elsewhere \cite{Joseph2017}. Briefly, agar powder was dissolved in de-ionised water and the solution was heated up to 100$^o$C inside a microwave oven. Pre-warmed intralipid and nigrosin dye solutions were added into the cooling agar solution to introduce scattering and absorption. A vortex mixer was used to homogenise the phantom solution and prevent separation of its components. The mixture was pipetted into petri dishes and the solutions were allowed to cool under room temperature to form phantoms having $\sim$2 mm thickness.

\subsection{Characterisation of optical absorption and scattering}
A double integrating sphere (DIS) system\cite{Pickering1993} was used to determine the reduced scattering and absorption coefficients of the phantoms. The DIS set-up uses two highly reflective spheres (Labsphere), with Lambertian surfaces made from polytetrafluoroethylene (PTFE), to capture and quantify the amount of transmitted and reflected light. A broadband tungsten halogen lamp (AvaLight-HAL-MINI, Avantes) and two fiber-optic spectrometers (AvaSpec-ULS2048, Avantes) served as the light source and detectors respectively. The inverse-adding doubling algorithm (IAD) was used to compute the optical properties of the phantoms from the reflectance, transmittance and reference measurements \cite{Prahl1993}.

\subsection{Characterisation of stretching deformation and associated birefringence}
Deformation contours across the birefringent phantoms mounted in the stretcher were assessed by drawing a grid pattern onto the polyacrylamide gel using a marker pen. Qualitative observation of the contours after stretching revealed a macroscopic non-linearity in the deformation, however, over the small field of view of the endoscope the stretch was approximately linear. Due to the constant volume of the material, the change in the thickness upon stretching must be accounted for when plotting retardance, by considering the change in area between the two clamps and dividing the measured retardance by the inferred width.

\pagebreak
\noindent
\Large\textbf{Supplementary Note 4:}\\
\Large\textbf{Calculation of entropy and expectation in endoscopic images}\\
\normalsize
\setcounter{section}{4}
\setcounter{subsection}{0}
\subsection{Distribution fitting}
To identify lesions in phase and polarisation images, two image filters/operators are used (see main manuscript): spatial entropy, $H()$, and spatial expectation, $E()$.  These can both determined under a unified framework that estimates the probability distribution of spatial points and then provides the differential entropy of the inferred distribution as $H()$ and the mean as $E()$.  This framework is also easily extended to joint distributions, e.g. for co-occurrence metrics or for multivariate entropy.

The distribution fitting process is detailed in Supplementary Figure \ref{fig:entropyCalcAlg}.  The first step is to defined a spatial region (a set of $M$ pixels, $\mathbf{S}$) to which a distribution will be fitted. It must be decided whether or not the pixel values $\mathbf{S}$ are assumed to be statistically independent.  In the most general case, the $M$ pixels of $\mathbf{S}$, will have a joint distribution density function $p(\mathbf{z}_1, ..., \mathbf{z}_m, ..., \mathbf{z}_M)$, where $\mathbf{z}_m$ is a $P$-dimensional vector containing some subset of available quantities (amplitude, phase, retardance etc.).  Under this model, all $M$ pixels together effectively represent a single sample from an $M\cdot P$-dimensional distribution.  It is not typically feasible to fit such a distribution based on only a single sample and so some assumptions must be made to make the problem tractable.

A widely used assumption is that points separated by some fixed spatial vector $(a,b)^T$ are drawn from a joint distribution, $p(\mathbf{z}(x,y), \mathbf{z}(x+a,y+b))$.  This only has $2\cdot P$ dimensions and there are now $M$ samples for fitting (one for each point in $\mathbf{S}$). This assumption is regularly used in grayscale image amplitude analysis where it is termed the \emph{gray-level co-occurrence matrx (GLCM)} \cite{Albregtsen2008}.  In this work we therefore refer to quantities generated using this assumption as being \emph{co-occurrence} metrics.

A further relaxation from the co-occurrence case is to assume that each pixel, $\mathbf{z}$ in the area $\mathbf{S}$ is drawn independently from a distribution, $p(\mathbf{z})$. This, too, offers $M$ points for fitting.  However, as discussed in the main manuscript, it is found that there is relatively little difference between these two cases for the samples measured and so the latter is typically preferred.

If the quantitites chosen for fitting are all circular (i.e. represent angles), then a \emph{generalisd multivariate von Mises distribution} is used for fitting \cite{Mardia2008}. If only $A$ is chosen a chi distribution is fitted ($A$ is modelled as the absolute value of a Gaussian) and if only $D$ is chosen, a truncated Gaussian on the range $[-1,1]$ is fitted.  If a mixture of circular and non-circular quantities is used it is difficult in general to define a multivariate distribution with appropriate marginal distributions for all quantities. One solution is to scale the domain of the non-circular quantities such that their expected probability mass lies close to 0 on the $(-\pi,\pi]$ range.  In this case, as we let $\kappa \rightarrow 0$, the von Mises distribution starts to approximate a truncated Gaussian distribution which can be fit to the non-circular quantities.

Fitting is performed by evaluating predetermined analytic expressions for maximum likelihood parameters of distributions based on measured quantities. Such expressions are known for Gaussian and multivariate von Mises distributions \cite{Mardia2008}.

\subsection{Computation of entropy}
\label{subsec:entropyCalc}
Fitting distributions to points in a spatial region gives estimates for the mean and variance of a quantity in that region.  Given these, the \emph{entropy} of the distribution can then be computed, which gives a measure of spatial heterogeneity that can be indicative of disordered microstructure. There are many different definitions of entropy (e.g. Shannon entropy, thermodynamic entropy), although on some level they are fundamentally analogous. The definition used here is the \emph{differential entropy}, which is a property of the fitted distributions.  This metric can be derived from the Kullback-Leibler divergence, which measures the similarity of probability distributions $P$ and $Q$ with density functions $p(x)$ and $q(x)$ (the 1-D case) respectively:

\begin{equation}
\label{eq:KLdiv}
D_{KL}(P||Q) = \int_{-\infty}^{\infty} p(x) \log \frac{p(x)}{q(x)} dx
\end{equation}

If $Q$ is a uniform distribution then $D_{KL}(P||Q)$ represents how similar to a uniform distribution $P$ is -- that is, how 'spread out' $P$ is.  This is precisely what we are seeking: a measure of heterogeneity.  Setting $Q$ to be a uniform distribution and introducing a scale factor of -1, Equation \ref{eq:KLdiv} becomes the differential entropy, $H$:

\begin{equation}
\label{eq:shannonEnt}
H = - \int_{-\infty}^{\infty} p(x) \log p(x) dx
\end{equation}

This can be extended to multivariate distributions simply by integrating over the additional variables.  However, computation of this integral can be computationally expensive.  Fortunately, for many distributions analytical expressions for the entropy in terms of distribution parameters are known. This is true of the multivariate Gaussian and univariate von Mises distributions \cite{Mardia2008}, shown in Supplementary Figure \ref{fig:entropyCalcAlg}.  Where analytical expressions are not known (e.g. multivariate von Mises) Monte-Carlo sampling methods are used for efficient integral evaluation \cite{MacKay2003}

By computing entropy for a sliding window over the entire raw image, we produce an \emph{entropy filtered} image of the desired quantity. Entropy filtering has been successfully applied in the field of amplitude-only texture analysis for classification of diseased tissues \cite{Yogesan1996,DeArruda2013,Ganeshan2007}.

\subsection{Comparison of entropy and expectation}
If instead of plotting entropy as described in the previous section, we plot the mean of the fitted distributions at each spatial point, this gives images of the spatial expectation, $E$, effectively an averaging filtering on the same length scale as entropy. As shown in Supplementary Figure \ref{fig:CNRall} it was found experimentally that phase entropy is a good indicator of early tumourigenesis, whereas for polarisation properties, expectation offers better performance.

In the case of retardance, $\eta$, this could be because diseased tissue has a higher \emph{average} retardance due to, for example, increasingly dense collagen networks associated with early malignant transformation. The spatial expectation would then represent this change better than entropy.  Further, the anisotropic nature of polarising molecules such as collagen may mean they are preferentially aligned along some axis.  Even if this alignment is weak, we might then expect to see an increase in the average value of the retardance axis orientation, $\theta_{\eta}$ or diattenuation axis orientation, $\theta_D$.  This is consistent with the results of Supplementary Figure \ref{fig:CNRall} that show more contrast in the spatial expectation of $\theta_{\eta}$ and $\theta_{\eta}$ quantities than in the entropy.

\pagebreak
\noindent
\Large\textbf{Supplementary Note 5:}\\
\Large\textbf{Limitations of current design and routes to real-time application}\\
\normalsize
\setcounter{section}{5}
\setcounter{subsection}{0}

\subsection{Experimental imaging acquisition speed}
Currently, the acquisition time for an amplitude, phase and polarisation image is 8.3 s. This is because 11 images are required: 7 for establishing the phase of the horizontally polarised image and 4 for determining the phase difference between the horizontally and vertically polarised images. To infer polarimetric properties of a tissue sample, 3 such images are recorded (see Online Methods and Supplementary Note \ref{subsec:bayesInf}) taking a total of 24.9 s.  For real time operation, this would need to be reduced substantially to $<$0.2 s.

The following modifications could be made to progress the operation towards real-time imaging. Using a camera with an automated high dynamic range imaging procedure, region of interest readout, and consequently higher frame rate could offer a >15 -fold increase in speed. Next, using an appropriately designed phase-mask in place of the current Fresnel lens, it would be possible to use only 2 images for phase retrieval and 3 for polarimetry, offering a further 2-fold speed up \cite{Horisaki2015,Zhang2007}.  Using only 2 illumination states instead of 3 offers a further reduction. Together these improvements would bring sample imaging time down to below the 0.2 s target. Further speed increases could be achieved through parallelisation of phase/polarisation retrieval by subdividing the camera sensor area, offering a trade-off with resolution \cite{Millerd2004}. At this point the switching speed of the liquid crystals in the SLM (60Hz frame rate) would become the limiting factor but this can be addressed by using high-speed (20kHz frame rate) digital micromirror devices (DMDs) \cite{Dremeau2015, Mitchell2016}.

\subsubsection{Transmission matrix recovery speed}
Although the fibre transmission matrix is typically stable for many months on an optical bench, in reality perturbations such as fibre bending would require periodic re-measurement of the transmission matrix. Our parallel spot characterisation system (described in Online Methods) already offers a speed-up of 12 over conventional spot or plane-wave scanning systems for recording fibre characterisation measurements. As opposed to methods used for MMFs, this approach scales to higher resolutions without requiring additional experimental measurements.

Currently, calibration measurements take 50.8 minutes but by making use of the optimisations discussed in the previous section this could be reduced to around 22 seconds. Re-characterisation may only need to be performed periodically at longer intervals, meaning it should not drastically affect total imaging time.  Further, parallelised compressed imaging techniques could be exploited with an array of high-speed photodiodes to move towards real-time characterisation if required \cite{Macfaden2016}.

\subsection{Image reconstruction speed}
To recover phase and polarisation information from the raw intensity images produced by the camera, a phase retrieval algorithm is used (see Online Methods). This currently takes around 67 s per single image (67 $\times$ 3 = 201 for full polarimetric imaging). The iterative algorithm used here (based on the Gerchberg-Saxton algorithm) requires around 200 iterations to converge. If we were to instead directly solve transport of intensity equations (TIE) this could be reduced to $<$0.1 s per single image \cite{Allen2001}.

For transmission matrix recovery, it is necessary to solve 293 phase retrieval problems (81$\times$3 characterisation $+$ 50 phase stabilisation), which currently takes 5.5 hours. Again by using a TIE approach this could be reduced to 29.5 s.  Because these problems are highly parallelisable, this could be sped up another $\sim$100-fold by using state-of-the-art processors (e.g. Intel Xeon server with 8 CPUs with 24 cores each).  Further parallelisation could be achieved by using dedicated graphics processing units (GPUs). The actual recovery of the matrix from these images requires solving 81$\times$12$\times$2=1944 inverse problems -- one per row of the desired $\mathbf{A}^{-1}$ (Supplementary Note \ref{subsec:tMatRecovery}).  This currently takes 9.5 hours.  However, by downsampling the camera images (e.g. sampling one pixel for each fibrelet) and increasing the degree of parallelisaiton from the current factor of 8 to $>$100, it should be possible to reduce this to $\sim$30 s. This could be further reduced by adopting solution methods that directly reconstruct input fields as weighted sums of the calibration fields rather than inverting the full transmission matrix \cite{Choi2014}.

Bayesian inference of sample polarimetric properties (described in Supplementary Note \ref{subsec:bayesInf}) is performed using the STAN package maximum likelihood optimiser, which requires 1.6 seconds to perform inference on a single pixel.  A full endoscope frame then takes 24 minutes to reconstruct.  Again, this problem is highly parallelisable and so could benefit from a 100-fold increase in speed on a state-of-the-art CPU. Further, by deriving analytical expressions for the joint maximum likelihood parameters of neighbouring pixels and using these instead of optimisation, inference time could be reduced to $<$0.2 s per image.

Computing the spatial entropy and expectation of an image, which involves fitting distributions as described in Supplementary Note \ref{subsec:entropyCalc}, currently takes 5 s per endoscope frame. Again, by utilising greater parallelisation this could be reduced to $<$0.1 s.

Finally, there are solution methods that do not require an explicit phase-retrieval step and can retrieve the desired transfer matrix based purely on intensity-only measurements. For the case where the matrix is known to comprise independent, complex Gaussian distributed entries a variational inference approach can be used \cite{Dremeau2015}. For the more generalised case in which the matrix is sparse, as it is here, algorithms such as \emph{generalised approximate message passing (GAMP)} may be applied directly to amplitude image sets \cite{Rangan2010, Schniter2015}. These could potentially reduce reconstruction time, as well as requiring minimal experimental measurements.

\subsection{Reflection mode operation}
For \emph{in vivo} clinical use, the endoscope must operate entirely in reflection mode rather than its current transmission mode.  In principle, the transmission mode illumination of samples can be simulated in reflection by using oblique back-scattering illumination \cite{Ford2012}. However, the imaging fibre will invariably experience random bending and temperature variations when in use and so the transmission matrix will need to be re-characterised without access to the distal facet.

One approach to reflection mode operation is to attach a calibration beacon to the distal facet, which has been demonstrated to work for creating a single focussed spot in confocal imaging but does not in the same form support wide-field imaging \cite{Farahi2013}. Speckle correlation imaging does not require any calibration of fibres, but relies on the assumption that images are of real, positive targets with finite support, which is not true in phase imaging where objects are generally complex \cite{Porat2016}. Another approach is to attach a known reflective plate to the distal facet and use back reflected light to make dynamic adjustments to a pre-recorded transmission matrix \cite{Gu2015}, however, this has not yet been demonstrated experimentally. Initial work using highly-spaced MCF has shown promise for two-photon imaging \cite{Warren2016}. Experimental investigation of these approaches is underway in our laboratories as a next step towards clinical translation.

\input{supRefs.bbl}
\bibliographystyle{naturemag}
\end{document}