Let $J\left(U_{0},\ldots,U_{N-1}\right)=\sum_{k=0}^{N}trace\left(X_{k}^{*}X_{k}\right)$,

where $X_{0}$ is given and $X_{k+1}=A_{k}X_{k}+B_{k}U_{k}$, for $k=0,\ldots,N-1$, where $X_{k}$ is $m\times n$ matrix, $A_{k}$ is $m\times m$ matrix, $B_{k}$ is $m\times p$ matrix and $U_{k}$ is 0-1 $p\times n$ matrix, for $k=0,\ldots,N-1$.

The problem is: Find $\max J\left(U_{0},\ldots,U_{N-1}\right)$ uder the constraints: $X_{k+1}=A_{k}X_{k}+B_{k}U_{k}$ and $U_{k}$ are 0-1 matrices, for $k=0,\ldots,N-1$. Also, find a sequence $U_{0},\ldots,U_{N-1}$ of 0-1 matrices, for which $\max J\left(U_{0},\ldots,U_{N-1}\right)$ is acheived.

Unfolding the dynamics:

\begin{equation}

\begin{aligned}

X_{1}&=A_{0}X_{0}+B_{0}U_{0}\\

X_{2}&=A_{1}X_{1}+B_{1}U_{1}\\

&\vdots\\

X_{N}&=A_{N-1}X_{N-1}+B_{N-1}U_{N-1},

\end{aligned}

\end{equation}

implies

\begin{equation}

\begin{aligned}

X_{1}&=A_{0}X_{0}+B_{0}U_{0}\\

X_{2}&=A_{1}A_{0}X_{0}+A_{1}B_{0}U_{0}+B_{1}U_{1}\\

&=A_{1:0}X_{0}+A_{1}B_{0}U_{0}+B_{1}U_{1}\\

&\vdots\\

X_{k}&=A_{k-1:0}X_{0}+A_{k-1:1}B_{0}U_{0}+A_{k-1:2}B_{1}U_{1}+\cdots+A_{k-1:k-1}B_{k-2}U_{k-2}+B_{k-1}U_{k-1}\\

&\vdots\\

X_{N}&=A_{N-1:0}X_{0}+A_{B-1:1}B_{0}U_{0}+A_{N-1:2}B_{1}U_{1}+\cdots+A_{N-1:N-1}B_{T-2}U_{N-2}+B_{N-1}U_{N-1},

\end{aligned}

\end{equation}

where $A_{i:j}=A_{i}A_{i+1}\cdots A_{j}$, for any $0\leq i\leq j\leq N-1$.

Writing the last in matrix form, it follows that

\begin{equation}

\begin{aligned}

&\left[\begin{array}{c}X_{1}\\X_{2}\\\vdots\\X_{k}\\\vdots\\X_{N}\end{array}\right]=\left[\begin{array}{c}A_{0:0}\\A_{1:0}\\\vdots\\A_{k-1:0}\\\vdots\\A_{N-1:0}\end{array}\right]X_{0}+\\

&+\left[\begin{array}{cccccccc}

B_{0}&0&0&0&0&0&\cdots&0\\

A_{1:1}B_{0}&B_{1}&0&0&0&0&\cdots&0\\

\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\

A_{k-1:1}B_{0}&A_{k-1:2}B_{1}&\cdots&A_{k-1:k-1}B_{k-2}&B_{k-1}&0&\cdots&0\\

\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\

A_{N-1:1}B_{0}&A_{N-1:2}B_{1}&\cdots&\cdots&\cdots&\cdots&A_{N-1:N-1}B_{N-2}&B_{N-1}

\end{array}\right]

\left[\begin{array}{c}U_{0}\\U_{1}\\\vdots\\U_{k-1}\\\vdots\\U_{N-1}\end{array}\right],

\end{aligned}

\end{equation}

which is written as

\begin{equation}

{\cal X}={\cal A}X_{0}+{\cal B}{\cal U}.

\end{equation}

One can now apply the Lagrange method (under the constraint that ${\cal U}$ is an $Np\times n$ $0-1$ matrix) in order to find an optimal solution. However, it was found that the Lagrange method reveals nothing but a brute-force search. Therefore, the constraint was relaxed to finding a solution on the boundary of the ball of matrices with center $\frac{1}{2}{\bf 1}_{Np\times n}$ and radius $\frac{1}{2}\sqrt{Npn}$ (where ${\bf 1}_{Np\times n}$ denotes the $Np\times n$ matrix of $1$'s). The Lagrange method has now revealed an optimal solution that can be approximated in polynomial time. The solution is given as

\begin{equation}

{\cal U}\left(\lambda\right)=\left(\lambda I_{Ns}+{\cal B}^{*}{\cal B}\right)^{-1}\left(\frac{1}{2}\lambda{\bf 1}_{Ns\times n}-{\cal B}^{*}{\cal A}X_{0}\right),

\end{equation}

where $\lambda >0$, and an approximation can be found by applying Newton's method to find $\lambda>0$ for which $\left\|{\cal U}\left(\lambda\right)-\frac{1}{2}\lambda{\bf 1}_{Ns\times n}\right\|_{F}=\frac{1}{2}\sqrt{Npn}$.

My question is: How can Dynamic Programming be utilized to reduce the complexity further?