Control Theory \(101-\varepsilon\)

Control System

• Given: system governed by diff-eq \[x'(t)=\mathsf{Model}(t,x(t), c(t))\] and observations \[y(t)=\mathsf{Obs}(t,x(t),c(t))\]
• Goal I: choose \(c(t)\) to "control" \(x(t)\) to desired states, at least as observed by \(y(t)\).
• Goal II: estimate \(x(t)\) from \(y(t)\).

Linear Time Independent (LTI) CS

\(\exists\) matrices \(A,B,C,D\) where \[x'(t)=\mathsf{Model}(t,x,c)=Ax(t)+Bc(t)\] \[y(t)=\mathsf{Obs}(t,x,c) = Cx(t)+Dc(t)\]

stabilizable if \(\exists\) matrix \(F\) \(c(t)=Fx(t)\).

(LTI) Control System

Assume: \(\exists\) matrices \(A,B,C,D\) where \[x'(t)=\mathsf{Model}(t,x,c)=Ax(t)+Bc(t)\] \[y(t)=\mathsf{Obs}(t,x,c) = Cx(t)+Dc(t)\]

stabilizable if \(\exists\) matrix \(F\) \(c(t)=Fx(t)\)

(LTI) Control System

LTI Control Systems

SOLVE FOR \(X\) IN \[HX-XA=-FC\]

Equation generalizes in shape hundreds of ways.

Really interested?

• Peter Benner's Control Theory Basics, in Handbook of Linear Algebra, Chapman & Hall, 2007

Solving \(XA-BY=C\)

Simplifying assumptions

• matrices are all approximately \(d\times d\)
• field operations are exact (numerical stability not considered here)

Case: \(B\) invertible

\[\begin{aligned} XA-BY &= C\\ Y & = B^{-1}(C+XA)\end{aligned}\]

So every \(X\) has a \(Y\).

Compute \(Y\) from \(X\) in \(O(d^3)\)-ops.

Case: \(B\) left-invertible

(tall-skinny \(B\) full rank)

\[\begin{aligned} XA-BY &= C\\ Y & = B^{\#}(C+XA)\end{aligned}\]

Again every \(X\) has a \(Y\).

Construct \(Y\) from given \(X\) in \(O(d^3)\)-ops.

Solving \(XA-BY=C\)

when \(X=Y\)?

Case: \(B\) invertible

\[\begin{aligned} XA-BX &= C & X & = B^{-1}(C+XA)\end{aligned}\]

...no clear solution.

Direct attack...

\[XA-BX = C\]

• System of linear equations: \[\begin{aligned}\sum_{\ell} X_{i\ell}A_{\ell k}-\sum_{m}B_{im}X_{mk}& =C_{ij}\\ (A^t\otimes I_d-I_d\otimes B)vec(X)& =vec(C)\end{aligned}\]
• \(d^2\) equations in \(d^2\) unknowns...\(O(d^6)\)-ops

d=1000 makes for exaflops of work!

100 Mega Watts @ 10 gigflops/Watt.

vs. 0.1 Watt for \(O(d^3)\) version.

Sylvester's Theorem: If \(A,B\) are square and have a common eigen value then \[(\exists ! X)(XA-BX=C).\]

By dimensions expect at most one solution, i.e. \(d^2\) equations in \(d^2\) unknowns.

Bartels-Stewart

Ways to get Spec(A)

Iteratively apply \(A\) to vector \(v\) till \(\exists\alpha_i\)\[\begin{aligned} A^{k+1} v& =\alpha_0 A^0v+\cdots +\alpha_k A^k v \\ ann_A(x) & = x^{k+1}-a_{k-1}x^{k-1}-\cdots -a_1 x^1 -a_0\end{aligned}\] Choices in how to (a) discover relation (b) apply relation.

• (Algebraic Lens) Krylov, Lanczos, Coppersmith, Weidemann,
• (Geometric Lens) QR, SVD,
• (Practical Lens) software has a button I can push.

Bartels-Stewart Theorem. There is an \(O(d^3)\)-time algorithm to solve \(XA-BX=C\).

Find Schur Decomposition

\[\begin{aligned} F^*AF & = \begin{bmatrix} I & * \\ 0 & \tilde{A}\end{bmatrix} & EBE^* & = \begin{bmatrix} I & * \\ 0 & \tilde{B}\end{bmatrix}\end{aligned}\]

Solve smaller system \[\tilde{X}\tilde{A}-\tilde{B}\tilde{X}=\tilde{C}\] and pull back through

\[EXAF+EBYF=ECF\]

Why did \(X=Y\) make solving so hard?

Why so few solutions?

\(X=Y\) shifts from one Sylvester system to a simultaneous system...

\[\begin{aligned} XA_1 -B_1 Y & = C_1\\ XA_2-B_2 Y & = C_2\end{aligned}\] where \(A_2=I_d, B_2=I_d,C_2=0\).

Same happens if \(X=X^t, \bar{X}^t\)

Real problem is to solve:

\[(\forall i)(XA_i - B_i Y=C_i)\]

\[(\forall i)(XA_i - B_i Y=C_i)\] in algebra

• Roth Removal Rule/Block-diagonalizing/Module direct sums \[\begin{bmatrix} I & -X \\ 0 & I \end{bmatrix}\begin{bmatrix}  B_i & C_i\\ 0 & A_i \end{bmatrix}\begin{bmatrix} I & X\\ 0 & I\end{bmatrix}=\begin{bmatrix} B_i & 0\\ 0 & A_i\end{bmatrix}\]
• Endomorphisms of modules \[End(_R M)=\{ X\mid (\forall i)(XR_i=R_iX)\}\]
• Adjoint algebras \[\{(X,Y)\mid (\forall i)(XA_i=B_iY)\}\]
• Schneider conjecture: \(O(d^6)\) is the complexity

\[\sum_{\ell} X_{i\ell} A_{\ell jk} -\sum_{m}B_{im}Y_{mjk} =C_{ijk}\]

• \(d^2 e\) equations in \(2d^2\) variables, highly over determined (expect few to now solutions)
• Not very sparse: 1/d nonzeros; 
• Random \(d^2\times d^2\) subsystem is basically random dense, so \(O(d^6)\)-ops to solve

(Technically this is the matrix for a transposed version.)

Tempted to echelonize?

It Shave an order of magnitude, but fills in quick.

And yet a ton of symmetry...exploit this!

Solving

\[(\forall i)(XA_i - B_i Y=C_i)\]

In \(O(d^3)\) ops away from edge cases

It is a tensor problem

• \([A]=[A_1,\ldots,A_c]\in K^s\otimes K^b\otimes K^c\)
• \([B]=[B_1,\ldots,B_c]\in K^a\otimes K^t\otimes K^c\)
• \([C]=[C_1,\ldots,C_c]\in K^a\otimes K^b\otimes K^c\)
• Now under tensor contraction (generalized matrix multiplication) solve one tensor equation \[X [A]-[B]Y=[C]\]

Solution Scheme Repeat Bartels-Stewart

• Reduce dimension via operators \(E,F\) \[\begin{aligned} EX [A]F+E[B]YF & =E[C]F \\ \tilde{X}[\tilde{A}]+[\tilde{B}]\tilde{Y} & = \tilde{C}\end{aligned}\]
• Reconstruct solution from lower-dimensional solution

What is different?

• \(E,F\) are 4-tensors, denote \([[E]]\) (was 2-tensors) \[\begin{aligned} [[E]]X [A][[F]]+[[E]][B]Y[[F]] & =[[E]][C][[F]] \\ \tilde{X}[\tilde{A}]+[\tilde{B}]\tilde{Y} & = \tilde{C}\end{aligned}\]
• These are tensor contractions and it is quite incorrect to write as binary products, e.g. [A][[F]], we need new notation.

2-tensor/Matrix

• for the dot-product \[\sum_{\ell} X_{i\ell} A_{\ell jk}\]
• more generally can replace dot-product with any bilinear form \(\langle|\rangle:U\times U\rightarrowtail K\) 

\([[E]]X[A][[F]]\)

\(-[[E]][B]Y[[F]]\)

\(=[[E]][C][[F]]\)

\[\begin{bmatrix} (A^1)^{-1} & \\ -(A^1)^{-1}A^2& I & & \\ \vdots & & \ddots & \\ -(A^1)^{-1} A^c & & & I\end{bmatrix}\begin{bmatrix} A^1\\ A^2 \\ \vdots \\ A^c\end{bmatrix}=\begin{bmatrix}A^1\\ 0\\ \vdots\\ 0\end{bmatrix}\]

Schur Decomposition in Matrix blowup

Schur Decomposition Tensorially

In tensor form...use face to clear the drawer (face elimination)

\[\begin{bmatrix} A_1^{-1} & \\ -A_1^{-1}A_2& I & & \\ \vdots & & \ddots & \\ -A_1^{-1} A_c & & & I\end{bmatrix}\begin{bmatrix} A_1\\ A_2 \\ \vdots \\ A_c\end{bmatrix}=\begin{bmatrix}A_1\\ 0\\ \vdots\\ 0\end{bmatrix}\]

In tensor form...use face to clear the drawer (face elimination)

In tensor form...use face to clear the drawer (face elimination)

=

\([[E]]X[A][[F]]\)

\(-[[E]][B]Y[[F]]\)

\(=[[E]][C][[F]]\)

(As [[E]] slides past [[F]] we can face echelonize [B])

(Face echelonize [A] leaves a tiny system to solve)