Consider the finite-dimensional continuous-time LTI realization

\[ \dot{\mathbf{x} }(t)=\mathbf{A}\mathbf{x}(t)+\mathbf{B}\mathbf{u}(t), \qquad \mathbf{y}(t)=\mathbf{C}\mathbf{x}(t)+\mathbf{D}\mathbf{u}(t). \]

The controllability matrix, observability matrix, reachable subspace, and unobservable subspace are

\[ \mathcal{C}= \begin{bmatrix}\mathbf{B} & \mathbf{A}\mathbf{B} & \cdots & \mathbf{A}^{n-1}\mathbf{B}\end{bmatrix},\qquad \mathcal{O}= \begin{bmatrix}\mathbf{C}\\ \mathbf{C}\mathbf{A}\\ \vdots\\ \mathbf{C}\mathbf{A}^{n-1}\end{bmatrix}, \]

\[ \mathscr{R}=\operatorname{im}\mathcal{C},\qquad \mathscr{N}=\ker\mathcal{O}. \]

The key invariance relations are

\[ \mathbf{A}\mathscr{R}\subseteq\mathscr{R},\qquad \mathbf{A}\mathscr{N}\subseteq\mathscr{N},\qquad \operatorname{im}\mathbf{B}\subseteq\mathscr{R},\qquad \mathscr{N}\subseteq\ker\mathbf{C}. \]

Kalman decomposition is the coordinate system in which these four inclusions become visible as zeros in the transformed matrices.

The state space is split into four parts:

1. Controllable-observable: reachable from input and visible at output.
2. Controllable-unobservable: reachable from input but hidden from output.
3. Uncontrollable-observable: not input-reachable but visible from initial conditions.
4. Uncontrollable-unobservable: neither input-reachable nor output-visible.

Choose subspaces satisfying

\[ \mathscr{X}_2=\mathscr{R}\cap\mathscr{N},\qquad \mathscr{R}=\mathscr{X}_1\oplus\mathscr{X}_2,\qquad \mathscr{N}=\mathscr{X}_2\oplus\mathscr{X}_4, \]

\[ \mathbb{R}^n = \mathscr{X}_1\oplus\mathscr{X}_2\oplus \mathscr{X}_3\oplus\mathscr{X}_4. \]

flowchart TD
  S["State space"] --> R["Reachable \nsubspace R"]
  S --> N["Unobservable \nsubspace N"]
  R --> X1["X1: controllable \nobservable"]
  R --> X2["X2: controllable \nunobservable"]
  N --> X2
  N --> X4["X4: uncontrollable \nunobservable"]
  S --> X3["X3: uncontrollable \nobservable"]
  X1 --> G["Contributes to \ntransfer matrix"]
  X2 --> H1["Input-reachable \nbut hidden"]
  X3 --> H2["Visible from initial \ncondition only"]
  X4 --> H3["Hidden autonomous \ninternal dynamics"]
        

Let \( \mathbf{T}=[\mathbf{T}_1\;\mathbf{T}_2\;\mathbf{T}_3\;\mathbf{T}_4] \), where the columns of \( \mathbf{T}_i \) span \( \mathscr{X}_i \), and define \( \mathbf{x}=\mathbf{T}\mathbf{z} \). Then

\[ \bar{\mathbf{A} }=\mathbf{T}^{-1}\mathbf{A}\mathbf{T},\qquad \bar{\mathbf{B} }=\mathbf{T}^{-1}\mathbf{B},\qquad \bar{\mathbf{C} }=\mathbf{C}\mathbf{T}. \]

With block order \( \mathbf{z}=[\mathbf{z}_1^T,\mathbf{z}_2^T,\mathbf{z}_3^T,\mathbf{z}_4^T]^T \), the Kalman decomposition has the structural form

\[ \bar{\mathbf{A} }= \begin{bmatrix} \mathbf{A}_{11} & \mathbf{0} & \mathbf{A}_{13} & \mathbf{0}\\ \mathbf{A}_{21} & \mathbf{A}_{22} & \mathbf{A}_{23} & \mathbf{A}_{24}\\ \mathbf{0} & \mathbf{0} & \mathbf{A}_{33} & \mathbf{0}\\ \mathbf{0} & \mathbf{0} & \mathbf{A}_{43} & \mathbf{A}_{44} \end{bmatrix},\qquad \bar{\mathbf{B} }= \begin{bmatrix}\mathbf{B}_1\\\mathbf{B}_2\\\mathbf{0}\\\mathbf{0}\end{bmatrix},\qquad \bar{\mathbf{C} }= \begin{bmatrix}\mathbf{C}_1&\mathbf{0}&\mathbf{C}_3&\mathbf{0}\end{bmatrix}. \]

The zero blocks are not arbitrary. They are forced by \( \mathbf{A}\mathscr{R}\subseteq\mathscr{R} \), \( \mathbf{A}\mathscr{N}\subseteq\mathscr{N} \), \( \operatorname{im}\mathbf{B}\subseteq\mathscr{R} \), and \( \mathscr{N}\subseteq\ker\mathbf{C} \).

Since \( \mathscr{R}=\mathscr{X}_1\oplus\mathscr{X}_2 \), the first and second block columns of \( \bar{\mathbf{A} } \) must map back into the first two block rows. Thus the third and fourth block rows in columns 1 and 2 are zero. Since \( \mathscr{N}=\mathscr{X}_2\oplus\mathscr{X}_4 \), the second and fourth block columns must map back into the second and fourth block rows. Thus the first and third block rows in columns 2 and 4 are zero. Combining these restrictions gives the matrix in Section 3.

The transformed transfer matrix is

\[ \mathbf{G}(s)= \bar{\mathbf{C} }(s\mathbf{I}-\bar{\mathbf{A} })^{-1}\bar{\mathbf{B} }+\mathbf{D}. \]

Under zero initial condition, only the controllable-observable block contributes:

\[ \mathbf{G}(s)= \mathbf{C}_1(s\mathbf{I}-\mathbf{A}_{11})^{-1}\mathbf{B}_1+\mathbf{D}. \]

Therefore the dimension of \( \mathbf{A}_{11} \) is the dimension of the minimal part of the realization.

Let \( r=\dim\mathscr{R} \), \( q=\dim\mathscr{N} \), and \( h=\dim(\mathscr{R}\cap\mathscr{N}) \). Then

\[ n_2=h,\qquad n_1=r-h,\qquad n_4=q-h,\qquad n_3=n-r-q+h. \]

These identities are useful for checking a numerical decomposition:

\[ r=\operatorname{rank}\mathcal{C},\qquad q=n-\operatorname{rank}\mathcal{O}. \]

A practical construction is:

1. Compute a basis of \( \mathscr{R}=\operatorname{im}\mathcal{C} \).
2. Compute a basis of \( \mathscr{N}=\ker\mathcal{O} \).
3. Compute \( \mathscr{R}\cap\mathscr{N} \) and call it \( \mathscr{X}_2 \).
4. Complete \( \mathscr{X}_2 \) to \( \mathscr{R} \) to obtain \( \mathscr{X}_1 \).
5. Complete \( \mathscr{X}_2 \) to \( \mathscr{N} \) to obtain \( \mathscr{X}_4 \).
6. Complete \( \mathscr{X}_1\oplus\mathscr{X}_2\oplus\mathscr{X}_4 \) to the full state space to obtain \( \mathscr{X}_3 \).

flowchart TD
  A["Compute Ctrb and basis R"] --> B["Compute Obsv and basis N"]
  B --> C["Intersection R cap N"]
  C --> D["X2: controllable unobservable"]
  D --> E["Complete to R: X1"]
  D --> F["Complete to N: X4"]
  E --> G["Complete full basis: X3"]
  F --> G
  G --> H["T = [X1 X2 X3 X4]"]
  H --> I["Abar, Bbar, Cbar"]
        

File: Chapter19_Lesson3.py

# Chapter19_Lesson3.py
# Kalman Decomposition: block structure of A, B, C
# Requires: numpy, scipy
#
# This script constructs a Kalman decomposition basis for
# x_dot = A x + B u, y = C x + D u.
#
# Block order used:
#   z = [ z_co ; z_cuo ; z_uco ; z_ucuo ]
# where
#   co   = controllable and observable part,
#   cuo  = controllable and unobservable part,
#   uco  = uncontrollable and observable part,
#   ucuo = uncontrollable and unobservable part.

import numpy as np
from scipy.linalg import svd, null_space

np.set_printoptions(precision=4, suppress=True)


def matrix_power_sequence(A, k):
    Ak = np.eye(A.shape[0])
    powers = []
    for _ in range(k):
        powers.append(Ak.copy())
        Ak = A @ Ak
    return powers


def controllability_matrix(A, B):
    n = A.shape[0]
    return np.hstack([Ak @ B for Ak in matrix_power_sequence(A, n)])


def observability_matrix(A, C):
    n = A.shape[0]
    return np.vstack([C @ Ak for Ak in matrix_power_sequence(A, n)])


def orth_basis(M, tol=1e-10):
    """Column-space orthonormal basis of M."""
    if M.size == 0:
        return np.zeros((M.shape[0], 0))
    U, s, _ = svd(M, full_matrices=False)
    r = np.sum(s > tol)
    return U[:, :r]


def rank(M, tol=1e-10):
    if M.size == 0:
        return 0
    return int(np.sum(svd(M, compute_uv=False) > tol))


def intersection_basis(U, V, tol=1e-10):
    """
    Basis for range(U) intersection range(V), assuming U,V have independent columns.
    We solve U a = V b, i.e. [U -V] [a;b] = 0.
    """
    if U.shape[1] == 0 or V.shape[1] == 0:
        return np.zeros((U.shape[0], 0))
    K = null_space(np.hstack([U, -V]), rcond=tol)
    if K.shape[1] == 0:
        return np.zeros((U.shape[0], 0))
    alpha = K[: U.shape[1], :]
    W = U @ alpha
    return orth_basis(W, tol)


def independent_append(current, candidates, target_dim=None, tol=1e-10):
    """
    Append independent columns from candidates to current.
    If target_dim is given, stop when that dimension is reached.
    """
    cols = []
    if current is not None and current.shape[1] > 0:
        cols = [current[:, j:j+1] for j in range(current.shape[1])]

    M = np.hstack(cols) if cols else np.zeros((candidates.shape[0], 0))
    current_rank = rank(M, tol)

    for j in range(candidates.shape[1]):
        c = candidates[:, j:j+1]
        trial = np.hstack([M, c])
        if rank(trial, tol) > current_rank:
            cols.append(c)
            M = trial
            current_rank += 1
            if target_dim is not None and current_rank >= target_dim:
                break

    return np.hstack(cols) if cols else np.zeros((candidates.shape[0], 0))


def complement_inside(container_basis, sub_basis, tol=1e-10):
    """Return columns that complete sub_basis to container_basis."""
    completed = independent_append(sub_basis, container_basis, target_dim=container_basis.shape[1], tol=tol)
    return completed[:, sub_basis.shape[1]:]


def kalman_decomposition(A, B, C, tol=1e-9):
    n = A.shape[0]

    R = orth_basis(controllability_matrix(A, B), tol)       # reachable subspace
    N = orth_basis(null_space(observability_matrix(A, C), rcond=tol), tol)  # unobservable subspace
    RN = intersection_basis(R, N, tol)                     # controllable-unobservable

    # V2: controllable and unobservable
    V2 = RN

    # V1 completes V2 to the controllable subspace: controllable-observable quotient
    V1 = complement_inside(R, V2, tol)

    # V4 completes V2 to the unobservable subspace: uncontrollable-unobservable quotient
    V4 = complement_inside(N, V2, tol)

    # V3 completes V1,V2,V4 to the full state space: uncontrollable-observable quotient
    I = np.eye(n)
    V124 = np.hstack([V1, V2, V4]) if (V1.shape[1]+V2.shape[1]+V4.shape[1]) else np.zeros((n,0))
    completed = independent_append(V124, I, target_dim=n, tol=tol)
    V3 = completed[:, V124.shape[1]:]

    # Kalman order: [co, c-unobs, unctrl-obs, unctrl-unobs]
    T = np.hstack([V1, V2, V3, V4])
    Tinv = np.linalg.inv(T)

    Abar = Tinv @ A @ T
    Bbar = Tinv @ B
    Cbar = C @ T

    dims = {
        "co": V1.shape[1],
        "c_unobs": V2.shape[1],
        "unctrl_obs": V3.shape[1],
        "unctrl_unobs": V4.shape[1],
        "rank_R": R.shape[1],
        "dim_N": N.shape[1],
    }
    return T, Abar, Bbar, Cbar, dims


def block_slices(dims):
    sizes = [dims["co"], dims["c_unobs"], dims["unctrl_obs"], dims["unctrl_unobs"]]
    idx = np.cumsum([0] + sizes)
    return [slice(idx[i], idx[i+1]) for i in range(4)]


def print_blocks(Abar, Bbar, Cbar, dims):
    labels = ["co", "c_unobs", "unctrl_obs", "unctrl_unobs"]
    slices = block_slices(dims)
    print("Block dimensions:", dims)
    print("\nAbar =")
    print(Abar)
    print("\nBbar =")
    print(Bbar)
    print("\nCbar =")
    print(Cbar)

    print("\nNorms of structurally zero A blocks:")
    zero_positions = [(0, 1), (0, 3), (2, 0), (2, 1), (2, 3), (3, 0), (3, 1)]
    for i, j in zero_positions:
        blk = Abar[slices[i], slices[j]]
        val = 0.0 if blk.size == 0 else np.linalg.norm(blk)
        print(f"  ||A_{labels[i]},{labels[j]}|| = {val:.2e}")

    print("\nNorms of structurally zero B/C blocks:")
    for i in [2, 3]:
        blk = Bbar[slices[i], :]
        print(f"  ||B_{labels[i]}|| = {np.linalg.norm(blk):.2e}")
    for j in [1, 3]:
        blk = Cbar[:, slices[j]]
        print(f"  ||C_{labels[j]}|| = {np.linalg.norm(blk):.2e}")


def demo():
    # Construct a system that already has the desired Kalman zero pattern.
    A_k = np.array([
        [-1.0,  0.0,  0.7,  0.0],
        [ 0.2, -2.0,  0.3, -0.4],
        [ 0.0,  0.0, -3.0,  0.0],
        [ 0.0,  0.0,  0.6, -4.0],
    ])
    B_k = np.array([[1.0],
                    [0.5],
                    [0.0],
                    [0.0]])
    C_k = np.array([[2.0, 0.0, -1.0, 0.0]])

    # Hide the block structure with a nonsingular physical-coordinate transform.
    T_phys = np.array([
        [1.0,  0.2,  0.1,  0.0],
        [0.0,  1.0, -0.3,  0.4],
        [0.2,  0.0,  1.0,  0.1],
        [0.1, -0.2,  0.0,  1.0],
    ])
    A = T_phys @ A_k @ np.linalg.inv(T_phys)
    B = T_phys @ B_k
    C = C_k @ np.linalg.inv(T_phys)

    _, Abar, Bbar, Cbar, dims = kalman_decomposition(A, B, C)
    print_blocks(Abar, Bbar, Cbar, dims)

    print("\nTransfer-function invariant check:")
    print("Only the controllable-observable block contributes to C(sI-A)^(-1)B + D.")


if __name__ == "__main__":
    demo()

File: Chapter19_Lesson3.cpp

// Chapter19_Lesson3.cpp
// Kalman Decomposition: block structure of A, B, C
// Requires Eigen: https://eigen.tuxfamily.org
//
// Compile example:
//   g++ -std=c++17 Chapter19_Lesson3.cpp -I /path/to/eigen -O2 -o Chapter19_Lesson3
//
// Block order:
//   [ controllable-observable,
//     controllable-unobservable,
//     uncontrollable-observable,
//     uncontrollable-unobservable ]

#include <Eigen/Dense>
#include <iostream>
#include <vector>
#include <cmath>

using Eigen::MatrixXd;
using Eigen::VectorXd;

double TOL = 1e-9;

int rankSVD(const MatrixXd& M, double tol = TOL) {
    if (M.size() == 0) return 0;
    Eigen::JacobiSVD<MatrixXd> svd(M);
    int r = 0;
    for (int i = 0; i < svd.singularValues().size(); ++i) {
        if (svd.singularValues()(i) > tol) r++;
    }
    return r;
}

MatrixXd nullSpace(const MatrixXd& M, double tol = TOL) {
    Eigen::JacobiSVD<MatrixXd> svd(M, Eigen::ComputeFullV);
    VectorXd s = svd.singularValues();
    int n = M.cols();
    int r = 0;
    for (int i = 0; i < s.size(); ++i) if (s(i) > tol) r++;
    int dim = n - r;
    if (dim <= 0) return MatrixXd(M.cols(), 0);
    return svd.matrixV().rightCols(dim);
}

MatrixXd columnSpaceBasis(const MatrixXd& M, double tol = TOL) {
    if (M.cols() == 0) return MatrixXd(M.rows(), 0);
    Eigen::JacobiSVD<MatrixXd> svd(M, Eigen::ComputeThinU);
    VectorXd s = svd.singularValues();
    int r = 0;
    for (int i = 0; i < s.size(); ++i) if (s(i) > tol) r++;
    if (r == 0) return MatrixXd(M.rows(), 0);
    return svd.matrixU().leftCols(r);
}

MatrixXd hstack(const std::vector<MatrixXd>& mats) {
    int rows = mats.empty() ? 0 : mats[0].rows();
    int cols = 0;
    for (const auto& M : mats) cols += M.cols();
    MatrixXd H(rows, cols);
    int c = 0;
    for (const auto& M : mats) {
        H.block(0, c, rows, M.cols()) = M;
        c += M.cols();
    }
    return H;
}

MatrixXd vstack(const std::vector<MatrixXd>& mats) {
    int cols = mats.empty() ? 0 : mats[0].cols();
    int rows = 0;
    for (const auto& M : mats) rows += M.rows();
    MatrixXd V(rows, cols);
    int r = 0;
    for (const auto& M : mats) {
        V.block(r, 0, M.rows(), cols) = M;
        r += M.rows();
    }
    return V;
}

MatrixXd controllabilityMatrix(const MatrixXd& A, const MatrixXd& B) {
    int n = A.rows();
    std::vector<MatrixXd> blocks;
    MatrixXd Ak = MatrixXd::Identity(n, n);
    for (int k = 0; k < n; ++k) {
        blocks.push_back(Ak * B);
        Ak = A * Ak;
    }
    return hstack(blocks);
}

MatrixXd observabilityMatrix(const MatrixXd& A, const MatrixXd& C) {
    int n = A.rows();
    std::vector<MatrixXd> blocks;
    MatrixXd Ak = MatrixXd::Identity(n, n);
    for (int k = 0; k < n; ++k) {
        blocks.push_back(C * Ak);
        Ak = A * Ak;
    }
    return vstack(blocks);
}

MatrixXd intersectionBasis(const MatrixXd& U, const MatrixXd& V, double tol = TOL) {
    if (U.cols() == 0 || V.cols() == 0) return MatrixXd(U.rows(), 0);
    MatrixXd M(U.rows(), U.cols() + V.cols());
    M << U, -V;
    MatrixXd K = nullSpace(M, tol);
    if (K.cols() == 0) return MatrixXd(U.rows(), 0);
    MatrixXd alpha = K.topRows(U.cols());
    return columnSpaceBasis(U * alpha, tol);
}

MatrixXd appendIndependent(const MatrixXd& current,
                           const MatrixXd& candidates,
                           int targetDim,
                           double tol = TOL) {
    MatrixXd M = current;
    int r = rankSVD(M, tol);
    for (int j = 0; j < candidates.cols(); ++j) {
        MatrixXd trial(M.rows(), M.cols() + 1);
        if (M.cols() > 0) trial.block(0, 0, M.rows(), M.cols()) = M;
        trial.col(M.cols()) = candidates.col(j);
        int rt = rankSVD(trial, tol);
        if (rt > r) {
            M = trial;
            r = rt;
            if (targetDim >= 0 && r >= targetDim) break;
        }
    }
    return M;
}

MatrixXd complementInside(const MatrixXd& container, const MatrixXd& sub, double tol = TOL) {
    MatrixXd completed = appendIndependent(sub, container, container.cols(), tol);
    if (completed.cols() <= sub.cols()) return MatrixXd(container.rows(), 0);
    return completed.rightCols(completed.cols() - sub.cols());
}

struct KalmanResult {
    MatrixXd T;
    MatrixXd Abar;
    MatrixXd Bbar;
    MatrixXd Cbar;
    std::vector<int> dims;
};

KalmanResult kalmanDecomposition(const MatrixXd& A, const MatrixXd& B, const MatrixXd& C) {
    int n = A.rows();

    MatrixXd R = columnSpaceBasis(controllabilityMatrix(A, B));
    MatrixXd N = columnSpaceBasis(nullSpace(observabilityMatrix(A, C)));
    MatrixXd V2 = intersectionBasis(R, N);
    MatrixXd V1 = complementInside(R, V2);
    MatrixXd V4 = complementInside(N, V2);

    MatrixXd V124 = hstack({V1, V2, V4});
    MatrixXd completed = appendIndependent(V124, MatrixXd::Identity(n, n), n);
    MatrixXd V3 = completed.rightCols(completed.cols() - V124.cols());

    MatrixXd T = hstack({V1, V2, V3, V4});
    MatrixXd Tinv = T.inverse();

    KalmanResult out;
    out.T = T;
    out.Abar = Tinv * A * T;
    out.Bbar = Tinv * B;
    out.Cbar = C * T;
    out.dims = {static_cast<int>(V1.cols()), static_cast<int>(V2.cols()),
                static_cast<int>(V3.cols()), static_cast<int>(V4.cols())};
    return out;
}

int main() {
    MatrixXd Ak(4,4);
    Ak << -1.0,  0.0,  0.7,  0.0,
           0.2, -2.0,  0.3, -0.4,
           0.0,  0.0, -3.0,  0.0,
           0.0,  0.0,  0.6, -4.0;

    MatrixXd Bk(4,1);
    Bk << 1.0, 0.5, 0.0, 0.0;

    MatrixXd Ck(1,4);
    Ck << 2.0, 0.0, -1.0, 0.0;

    MatrixXd Tphys(4,4);
    Tphys << 1.0,  0.2,  0.1,  0.0,
             0.0,  1.0, -0.3,  0.4,
             0.2,  0.0,  1.0,  0.1,
             0.1, -0.2,  0.0,  1.0;

    MatrixXd A = Tphys * Ak * Tphys.inverse();
    MatrixXd B = Tphys * Bk;
    MatrixXd C = Ck * Tphys.inverse();

    KalmanResult kr = kalmanDecomposition(A, B, C);

    std::cout << "Block dimensions [co, c_unobs, unctrl_obs, unctrl_unobs]: ";
    for (int d : kr.dims) std::cout << d << " ";
    std::cout << "\n\nAbar:\n" << kr.Abar << "\n\nBbar:\n" << kr.Bbar
              << "\n\nCbar:\n" << kr.Cbar << "\n";

    return 0;
}

File: Chapter19_Lesson3.java

// Chapter19_Lesson3.java
// Kalman Decomposition: block structure of A, B, C
//
// This pure-Java educational implementation uses rank tests and RREF-based
// bases. For production numerical work, use EJML, Apache Commons Math, or JBLAS.

import java.util.Arrays;

public class Chapter19_Lesson3 {
    static final double TOL = 1e-9;

    static class Result {
        double[][] T, Tinv, Abar, Bbar, Cbar;
        int[] dims;
    }

    static double[][] eye(int n) {
        double[][] I = new double[n][n];
        for (int i = 0; i < n; i++) I[i][i] = 1.0;
        return I;
    }

    static double[][] copy(double[][] A) {
        double[][] B = new double[A.length][A[0].length];
        for (int i = 0; i < A.length; i++) B[i] = Arrays.copyOf(A[i], A[i].length);
        return B;
    }

    static double[][] mmul(double[][] A, double[][] B) {
        int m = A.length, p = A[0].length, n = B[0].length;
        double[][] C = new double[m][n];
        for (int i = 0; i < m; i++)
            for (int k = 0; k < p; k++)
                for (int j = 0; j < n; j++)
                    C[i][j] += A[i][k] * B[k][j];
        return C;
    }

    static double[][] hstack(double[][]... mats) {
        int rows = mats[0].length;
        int cols = 0;
        for (double[][] M : mats) cols += M[0].length;
        double[][] H = new double[rows][cols];
        int c = 0;
        for (double[][] M : mats) {
            for (int i = 0; i < rows; i++)
                for (int j = 0; j < M[0].length; j++)
                    H[i][c + j] = M[i][j];
            c += M[0].length;
        }
        return H;
    }

    static double[][] vstack(double[][]... mats) {
        int cols = mats[0][0].length;
        int rows = 0;
        for (double[][] M : mats) rows += M.length;
        double[][] V = new double[rows][cols];
        int r = 0;
        for (double[][] M : mats) {
            for (int i = 0; i < M.length; i++)
                for (int j = 0; j < cols; j++)
                    V[r + i][j] = M[i][j];
            r += M.length;
        }
        return V;
    }

    static double[][] subCols(double[][] A, int start, int count) {
        double[][] S = new double[A.length][Math.max(count, 0)];
        for (int i = 0; i < A.length; i++)
            for (int j = 0; j < count; j++)
                S[i][j] = A[i][start + j];
        return S;
    }

    static double[][] negate(double[][] A) {
        double[][] B = copy(A);
        for (int i = 0; i < B.length; i++)
            for (int j = 0; j < B[0].length; j++)
                B[i][j] = -B[i][j];
        return B;
    }

    static double[][] transpose(double[][] A) {
        double[][] T = new double[A[0].length][A.length];
        for (int i = 0; i < A.length; i++)
            for (int j = 0; j < A[0].length; j++)
                T[j][i] = A[i][j];
        return T;
    }

    static double[][] rref(double[][] A) {
        double[][] R = copy(A);
        int rows = R.length, cols = R[0].length;
        int lead = 0;
        for (int r = 0; r < rows && lead < cols; r++) {
            int i = r;
            while (i < rows && Math.abs(R[i][lead]) < TOL) i++;
            if (i == rows) {
                lead++;
                r--;
                continue;
            }
            double[] tmp = R[r]; R[r] = R[i]; R[i] = tmp;

            double pivot = R[r][lead];
            for (int j = 0; j < cols; j++) R[r][j] /= pivot;

            for (int ii = 0; ii < rows; ii++) {
                if (ii == r) continue;
                double factor = R[ii][lead];
                for (int j = 0; j < cols; j++) R[ii][j] -= factor * R[r][j];
            }
            lead++;
        }
        return R;
    }

    static int rank(double[][] A) {
        double[][] R = rref(A);
        int rank = 0;
        for (int i = 0; i < R.length; i++) {
            boolean nonzero = false;
            for (int j = 0; j < R[0].length; j++)
                if (Math.abs(R[i][j]) > TOL) nonzero = true;
            if (nonzero) rank++;
        }
        return rank;
    }

    static double[][] independentColumns(double[][] candidates) {
        int m = candidates.length;
        double[][] current = new double[m][0];
        int r = 0;
        for (int j = 0; j < candidates[0].length; j++) {
            double[][] col = subCols(candidates, j, 1);
            double[][] trial = hstack(current, col);
            int rt = rank(trial);
            if (rt > r) {
                current = trial;
                r = rt;
            }
        }
        return current;
    }

    static double[][] appendIndependent(double[][] current, double[][] candidates, int targetDim) {
        double[][] M = current;
        int r = rank(M);
        for (int j = 0; j < candidates[0].length; j++) {
            double[][] col = subCols(candidates, j, 1);
            double[][] trial = hstack(M, col);
            int rt = rank(trial);
            if (rt > r) {
                M = trial;
                r = rt;
                if (targetDim >= 0 && r >= targetDim) break;
            }
        }
        return M;
    }

    static double[][] nullspace(double[][] A) {
        double[][] R = rref(A);
        int rows = R.length, cols = R[0].length;
        boolean[] pivot = new boolean[cols];
        int[] pivotRowForCol = new int[cols];
        Arrays.fill(pivotRowForCol, -1);

        for (int i = 0; i < rows; i++) {
            int pc = -1;
            for (int j = 0; j < cols; j++) {
                if (Math.abs(R[i][j]) > TOL) { pc = j; break; }
            }
            if (pc >= 0) {
                pivot[pc] = true;
                pivotRowForCol[pc] = i;
            }
        }

        int freeCount = 0;
        for (int j = 0; j < cols; j++) if (!pivot[j]) freeCount++;
        double[][] N = new double[cols][freeCount];
        int idx = 0;
        for (int f = 0; f < cols; f++) {
            if (pivot[f]) continue;
            N[f][idx] = 1.0;
            for (int p = 0; p < cols; p++) {
                if (pivot[p]) {
                    int row = pivotRowForCol[p];
                    N[p][idx] = -R[row][f];
                }
            }
            idx++;
        }
        return N;
    }

    static double[][] inverse(double[][] A) {
        int n = A.length;
        double[][] Aug = new double[n][2*n];
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) Aug[i][j] = A[i][j];
            Aug[i][n+i] = 1.0;
        }
        double[][] R = rref(Aug);
        double[][] Inv = new double[n][n];
        for (int i = 0; i < n; i++)
            for (int j = 0; j < n; j++)
                Inv[i][j] = R[i][n+j];
        return Inv;
    }

    static double[][] controllabilityMatrix(double[][] A, double[][] B) {
        int n = A.length;
        double[][] Ak = eye(n);
        double[][] R = new double[n][0];
        for (int k = 0; k < n; k++) {
            R = hstack(R, mmul(Ak, B));
            Ak = mmul(A, Ak);
        }
        return R;
    }

    static double[][] observabilityMatrix(double[][] A, double[][] C) {
        int n = A.length;
        double[][] Ak = eye(n);
        double[][] O = new double[0][C[0].length];
        for (int k = 0; k < n; k++) {
            double[][] block = mmul(C, Ak);
            O = (O.length == 0) ? block : vstack(O, block);
            Ak = mmul(A, Ak);
        }
        return O;
    }

    static double[][] intersectionBasis(double[][] U, double[][] V) {
        if (U[0].length == 0 || V[0].length == 0) return new double[U.length][0];
        double[][] K = nullspace(hstack(U, negate(V)));
        if (K[0].length == 0) return new double[U.length][0];
        double[][] alpha = new double[U[0].length][K[0].length];
        for (int i = 0; i < U[0].length; i++)
            for (int j = 0; j < K[0].length; j++)
                alpha[i][j] = K[i][j];
        return independentColumns(mmul(U, alpha));
    }

    static double[][] complementInside(double[][] container, double[][] sub) {
        double[][] completed = appendIndependent(sub, container, container[0].length);
        return subCols(completed, sub[0].length, completed[0].length - sub[0].length);
    }

    static Result kalman(double[][] A, double[][] B, double[][] C) {
        int n = A.length;
        double[][] R = independentColumns(controllabilityMatrix(A, B));
        double[][] N = independentColumns(nullspace(observabilityMatrix(A, C)));
        double[][] V2 = intersectionBasis(R, N);
        double[][] V1 = complementInside(R, V2);
        double[][] V4 = complementInside(N, V2);
        double[][] V124 = hstack(V1, V2, V4);
        double[][] completed = appendIndependent(V124, eye(n), n);
        double[][] V3 = subCols(completed, V124[0].length, completed[0].length - V124[0].length);

        double[][] T = hstack(V1, V2, V3, V4);
        double[][] Tinv = inverse(T);

        Result out = new Result();
        out.T = T;
        out.Tinv = Tinv;
        out.Abar = mmul(mmul(Tinv, A), T);
        out.Bbar = mmul(Tinv, B);
        out.Cbar = mmul(C, T);
        out.dims = new int[]{V1[0].length, V2[0].length, V3[0].length, V4[0].length};
        return out;
    }

    static void printMatrix(String name, double[][] A) {
        System.out.println(name + ":");
        for (double[] row : A) {
            for (double v : row) System.out.printf("%10.5f ", v);
            System.out.println();
        }
    }

    public static void main(String[] args) {
        double[][] A = {
            {-1.0, 0.0, 0.7, 0.0},
            { 0.2,-2.0, 0.3,-0.4},
            { 0.0, 0.0,-3.0, 0.0},
            { 0.0, 0.0, 0.6,-4.0}
        };
        double[][] B = { {1.0},{0.5},{0.0},{0.0} };
        double[][] C = { {2.0,0.0,-1.0,0.0} };

        Result kr = kalman(A, B, C);
        System.out.println("Block dimensions [co, c_unobs, unctrl_obs, unctrl_unobs] = "
                           + Arrays.toString(kr.dims));
        printMatrix("Abar", kr.Abar);
        printMatrix("Bbar", kr.Bbar);
        printMatrix("Cbar", kr.Cbar);
    }
}

MATLAB users may also compare the result with ss, ctrb, obsv, and minreal when the Control System Toolbox is available. The script below keeps the decomposition explicit.

File: Chapter19_Lesson3.m

% Chapter19_Lesson3.m
% Kalman Decomposition: block structure of A, B, C
%
% Block order:
%   [ controllable-observable,
%     controllable-unobservable,
%     uncontrollable-observable,
%     uncontrollable-unobservable ]

clear; clc;

A = [-1.0  0.0  0.7  0.0;
      0.2 -2.0  0.3 -0.4;
      0.0  0.0 -3.0  0.0;
      0.0  0.0  0.6 -4.0];

B = [1.0; 0.5; 0.0; 0.0];
C = [2.0 0.0 -1.0 0.0];

[T,Abar,Bbar,Cbar,dims] = kalman_decomposition_lesson(A,B,C);

disp('Block dimensions [co, c_unobs, unctrl_obs, unctrl_unobs]:');
disp(dims);
disp('Abar ='); disp(Abar);
disp('Bbar ='); disp(Bbar);
disp('Cbar ='); disp(Cbar);

% Optional verification of the input-output transfer function:
% If Control System Toolbox is available:
% sys_full = ss(A,B,C,0);
% sys_minimal = minreal(sys_full);
% tf(sys_minimal)

function [T,Abar,Bbar,Cbar,dims] = kalman_decomposition_lesson(A,B,C)
    n = size(A,1);
    R = orth(controllability_matrix_lesson(A,B));
    N = orth(null(observability_matrix_lesson(A,C)));

    V2 = intersection_basis_lesson(R,N);  % controllable and unobservable
    V1 = complement_inside_lesson(R,V2);  % controllable-observable quotient
    V4 = complement_inside_lesson(N,V2);  % uncontrollable-unobservable quotient

    V124 = [V1 V2 V4];
    completed = append_independent_lesson(V124, eye(n), n);
    V3 = completed(:, size(V124,2)+1:end); % uncontrollable-observable quotient

    T = [V1 V2 V3 V4];

    Abar = T \ (A*T);
    Bbar = T \ B;
    Cbar = C*T;

    dims = [size(V1,2), size(V2,2), size(V3,2), size(V4,2)];
end

function R = controllability_matrix_lesson(A,B)
    n = size(A,1);
    R = [];
    Ak = eye(n);
    for k = 1:n
        R = [R Ak*B]; %#ok<AGROW>
        Ak = A*Ak;
    end
end

function O = observability_matrix_lesson(A,C)
    n = size(A,1);
    O = [];
    Ak = eye(n);
    for k = 1:n
        O = [O; C*Ak]; %#ok<AGROW>
        Ak = A*Ak;
    end
end

function W = intersection_basis_lesson(U,V)
    if isempty(U) || isempty(V)
        W = zeros(size(U,1),0);
        return;
    end
    K = null([U -V]);
    if isempty(K)
        W = zeros(size(U,1),0);
        return;
    end
    alpha = K(1:size(U,2),:);
    W = orth(U*alpha);
end

function extra = complement_inside_lesson(container, sub)
    completed = append_independent_lesson(sub, container, size(container,2));
    extra = completed(:, size(sub,2)+1:end);
end

function M = append_independent_lesson(current, candidates, targetDim)
    M = current;
    r = rank(M);
    for j = 1:size(candidates,2)
        trial = [M candidates(:,j)]; %#ok<AGROW>
        if rank(trial) > r
            M = trial;
            r = r + 1;
            if r >= targetDim
                break;
            end
        end
    end
end

File: Chapter19_Lesson3.nb

Notebook[{
Cell["Chapter19_Lesson3.nb", "Title"],
Cell["Kalman Decomposition: Block Structure of A, B, C", "Subtitle"],
Cell["Block order: {controllable-observable, controllable-unobservable, uncontrollable-observable, uncontrollable-unobservable}.", "Text"],
Cell[BoxData[
 RowBox[{
  RowBox[{"ClearAll", "[", "\"Global`*\"", "]"}], ";"}]], "Input"],
Cell[BoxData[
 RowBox[{
  RowBox[{"controllabilityMatrix", "[", 
   RowBox[{"A_", ",", "B_"}], "]"}], ":=", 
  RowBox[{"ArrayFlatten", "[", 
   RowBox[{"{", 
    RowBox[{"Table", "[", 
     RowBox[{
      RowBox[{"MatrixPower", "[", 
       RowBox[{"A", ",", "k"}], "]"}], ".", "B"}], ",", 
      RowBox[{"{", 
       RowBox[{"k", ",", "0", ",", 
        RowBox[{
         RowBox[{"Length", "[", "A", "]"}], "-", "1"}]}], "}"}]}], "]"}], 
    "}"}], "]"}]}]], "Input"],
Cell[BoxData[
 RowBox[{
  RowBox[{"observabilityMatrix", "[", 
   RowBox[{"A_", ",", "C_"}], "]"}], ":=", 
  RowBox[{"Join", "@@", 
   RowBox[{"Table", "[", 
    RowBox[{
     RowBox[{"C", ".", 
      RowBox[{"MatrixPower", "[", 
       RowBox[{"A", ",", "k"}], "]"}]}], ",", 
     RowBox[{"{", 
      RowBox[{"k", ",", "0", ",", 
       RowBox[{
        RowBox[{"Length", "[", "A", "]"}], "-", "1"}]}], "}"}]}], "]"}]}]}]], "Input"],
Cell[BoxData[
 RowBox[{
  RowBox[{"basisColumns", "[", "M_", "]"}], ":=", 
  RowBox[{"Transpose", "[", 
   RowBox[{"Orthogonalize", "[", 
    RowBox[{"Transpose", "[", "M", "]"}], "]"}], "]"}]}]], "Input"],
Cell[BoxData[
 RowBox[{
  RowBox[{"intersectionBasis", "[", 
   RowBox[{"U_", ",", "V_"}], "]"}], ":=", 
  RowBox[{"Module", "[", 
   RowBox[{
    RowBox[{"{", 
     RowBox[{"K", ",", "alpha", ",", "m"}], "}"}], ",", 
    RowBox[{
     RowBox[{"m", "=", 
      RowBox[{"Dimensions", "[", "U", "]"}]}], ";", 
     RowBox[{"If", "[", 
      RowBox[{
       RowBox[{
        RowBox[{"m", "[", 
         RowBox[{"[", "2", "]"}], "]"}], "==", "0"}], ",", 
       RowBox[{"Return", "[", 
        RowBox[{"ConstantArray", "[", 
         RowBox[{"0", ",", 
          RowBox[{"{", 
           RowBox[{
            RowBox[{"m", "[", 
             RowBox[{"[", "1", "]"}], "]"}], ",", "0"}], "}"}]}], "]"}], 
        "]"}]}], "]"}], ";", 
     RowBox[{"K", "=", 
      RowBox[{"NullSpace", "[", 
       RowBox[{"Transpose", "[", 
        RowBox[{"Join", "[", 
         RowBox[{
          RowBox[{"Transpose", "[", "U", "]"}], ",", 
          RowBox[{"Transpose", "[", 
           RowBox[{"-", "V"}], "]"}]}], "]"}], "]"}], "]"}]}], ";", 
     RowBox[{"If", "[", 
      RowBox[{
       RowBox[{"K", "==", 
        RowBox[{"{", "}"}]}], ",", 
       RowBox[{"ConstantArray", "[", 
        RowBox[{"0", ",", 
         RowBox[{"{", 
          RowBox[{
           RowBox[{"m", "[", 
            RowBox[{"[", "1", "]"}], "]"}], ",", "0"}], "}"}]}], "]"}], ",", 
       RowBox[{
        RowBox[{"alpha", "=", 
         RowBox[{"Transpose", "[", 
          RowBox[{"K", "[", 
           RowBox[{"[", 
            RowBox[{"All", ",", 
             RowBox[{"1", ";;", 
              RowBox[{"m", "[", 
               RowBox[{"[", "2", "]"}], "]"}]}]}], "]"}], "]"}], "]"}]}], ";", 
        RowBox[{"basisColumns", "[", 
         RowBox[{"U", ".", "alpha"}], "]"}]}]}], "]"}]}]}], "]"}]}]], "Input"],
Cell[BoxData[
 RowBox[{
  RowBox[{"appendIndependent", "[", 
   RowBox[{"current_", ",", "candidates_", ",", "target_"}], "]"}], ":=", 
  RowBox[{"Module", "[", 
   RowBox[{
    RowBox[{"{", 
     RowBox[{"M", "=", "current", ",", "r", ",", "trial"}], "}"}], ",", 
    RowBox[{
     RowBox[{"r", "=", 
      RowBox[{"MatrixRank", "[", "M", "]"}]}], ";", 
     RowBox[{"Do", "[", 
      RowBox[{
       RowBox[{
        RowBox[{"trial", "=", 
         RowBox[{"Join", "[", 
          RowBox[{"M", ",", 
           RowBox[{"candidates", "[", 
            RowBox[{"[", 
             RowBox[{"All", ",", 
              RowBox[{"{", "j", "}"}]}], "]"}], "]"}], ",", "2"}], "]"}]}], ";", 
        RowBox[{"If", "[", 
         RowBox[{
          RowBox[{
           RowBox[{"MatrixRank", "[", "trial", "]"}], ">", "r"}], ",", 
          RowBox[{
           RowBox[{"M", "=", "trial"}], ";", 
           RowBox[{"r", "++"}], ";", 
           RowBox[{"If", "[", 
            RowBox[{
             RowBox[{"r", ">=", "target"}], ",", 
             RowBox[{"Break", "[", "]"}]}], "]"}]}]}], "]"}]}], ",", 
       RowBox[{"{", 
        RowBox[{"j", ",", 
         RowBox[{"Dimensions", "[", "candidates", "]"}], "[", 
         RowBox[{"[", "2", "]"}], "]"}], "}"}]}], "]"}], ";", "M"}]}], "]"}]}]], "Input"],
Cell[BoxData[
 RowBox[{
  RowBox[{"A", "=", 
   RowBox[{"(", GridBox[{
      {"-1.0", "0.0", "0.7", "0.0"},
      {"0.2", "-2.0", "0.3", "-0.4"},
      {"0.0", "0.0", "-3.0", "0.0"},
      {"0.0", "0.0", "0.6", "-4.0"}
     }], ")"}]}], ";", 
  RowBox[{"B", "=", 
   RowBox[{"(", GridBox[{ {"1.0"}, {"0.5"}, {"0.0"}, {"0.0"} }], ")"}]}], ";", 
  RowBox[{"Cmat", "=", 
   RowBox[{"(", GridBox[{ {"2.0", "0.0", "-1.0", "0.0"} }], ")"}]}], ";"}]], "Input"],
Cell[BoxData[
 RowBox[{
  RowBox[{"R", "=", 
   RowBox[{"basisColumns", "[", 
    RowBox[{"controllabilityMatrix", "[", 
     RowBox[{"A", ",", "B"}], "]"}], "]"}]}], ";", 
  RowBox[{"Nsub", "=", 
   RowBox[{"Transpose", "[", 
    RowBox[{"NullSpace", "[", 
     RowBox[{"observabilityMatrix", "[", 
      RowBox[{"A", ",", "Cmat"}], "]"}], "]"}], "]"}]}], ";", 
  RowBox[{"V2", "=", 
   RowBox[{"intersectionBasis", "[", 
    RowBox[{"R", ",", "Nsub"}], "]"}]}], ";", 
  RowBox[{"V1", "=", 
   RowBox[{
    RowBox[{"appendIndependent", "[", 
     RowBox[{"V2", ",", "R", ",", 
      RowBox[{"Dimensions", "[", "R", "]"}], "[", 
      RowBox[{"[", "2", "]"}], "]"}], "]"}], "[", 
    RowBox[{"[", 
     RowBox[{"All", ",", 
      RowBox[{
       RowBox[{
        RowBox[{"Dimensions", "[", "V2", "]"}], "[", 
        RowBox[{"[", "2", "]"}], "]"}], "+", "1"}], ";;"}], "]"}], "]"}]}], 
  ";"}]], "Input"],
Cell[BoxData[
 RowBox[{
  RowBox[{"V4", "=", 
   RowBox[{
    RowBox[{"appendIndependent", "[", 
     RowBox[{"V2", ",", "Nsub", ",", 
      RowBox[{"Dimensions", "[", "Nsub", "]"}], "[", 
      RowBox[{"[", "2", "]"}], "]"}], "]"}], "[", 
    RowBox[{"[", 
     RowBox[{"All", ",", 
      RowBox[{
       RowBox[{
        RowBox[{"Dimensions", "[", "V2", "]"}], "[", 
        RowBox[{"[", "2", "]"}], "]"}], "+", "1"}], ";;"}], "]"}], "]"}]}], 
  ";", 
  RowBox[{"V124", "=", 
   RowBox[{"Join", "[", 
    RowBox[{"V1", ",", "V2", ",", "V4", ",", "2"}], "]"}]}], ";", 
  RowBox[{"completed", "=", 
   RowBox[{"appendIndependent", "[", 
    RowBox[{"V124", ",", 
     RowBox[{"IdentityMatrix", "[", 
      RowBox[{"Length", "[", "A", "]"}], "]"}], ",", 
     RowBox[{"Length", "[", "A", "]"}]}], "]"}]}], ";", 
  RowBox[{"V3", "=", 
   RowBox[{"completed", "[", 
    RowBox[{"[", 
     RowBox[{"All", ",", 
      RowBox[{
       RowBox[{
        RowBox[{"Dimensions", "[", "V124", "]"}], "[", 
        RowBox[{"[", "2", "]"}], "]"}], "+", "1"}], ";;"}], "]"}], "]"}]}], 
  ";"}]], "Input"],
Cell[BoxData[
 RowBox[{
  RowBox[{"T", "=", 
   RowBox[{"Join", "[", 
    RowBox[{"V1", ",", "V2", ",", "V3", ",", "V4", ",", "2"}], "]"}]}], ";", 
  RowBox[{"Abar", "=", 
   RowBox[{
    RowBox[{"Inverse", "[", "T", "]"}], ".", "A", ".", "T"}]}], ";", 
  RowBox[{"Bbar", "=", 
   RowBox[{
    RowBox[{"Inverse", "[", "T", "]"}], ".", "B"}]}], ";", 
  RowBox[{"Cbar", "=", 
   RowBox[{"Cmat", ".", "T"}]}], ";", 
  RowBox[{"{", 
   RowBox[{"MatrixForm", "[", "Abar", "]"}], ",", 
    RowBox[{"MatrixForm", "[", "Bbar", "]"}], ",", 
    RowBox[{"MatrixForm", "[", "Cbar", "]"}]}], "}"}]}]], "Input"]
},
WindowSize->{1200, 760},
WindowTitle->"Chapter19_Lesson3.nb"
]

Consider

\[ \mathbf{A}= \begin{bmatrix} -1 & 0 & 0.7 & 0\\ 0.2 & -2 & 0.3 & -0.4\\ 0 & 0 & -3 & 0\\ 0 & 0 & 0.6 & -4 \end{bmatrix},\quad \mathbf{B}=\begin{bmatrix}1\\0.5\\0\\0\end{bmatrix},\quad \mathbf{C}=\begin{bmatrix}2&0&-1&0\end{bmatrix}. \]

The transfer matrix is determined only by \( A_{11}=-1 \), \( B_1=1 \), and \( C_1=2 \):

\[ G(s)=2(s+1)^{-1}=\frac{2}{s+1}. \]

The remaining modes are internal nonminimal modes for the input-output map, although they may still affect initial-condition response.

Problem 1: Prove that \( \mathscr{R}=\operatorname{im}\mathcal{C} \) is \( \mathbf{A} \)-invariant.

Solution: Every reachable vector can be written as \( \mathbf{v}=\sum_{k=0}^{n-1}\mathbf{A}^k\mathbf{B}\boldsymbol{\alpha}_k \). Then

\[ \mathbf{A}\mathbf{v}= \sum_{k=0}^{n-1}\mathbf{A}^{k+1}\mathbf{B}\boldsymbol{\alpha}_k. \]

The term \( \mathbf{A}^n\mathbf{B} \) is a linear combination of lower powers by Cayley-Hamilton, so \( \mathbf{A}\mathbf{v}\in\mathscr{R} \).

Problem 2: Prove that \( \mathscr{N}=\ker\mathcal{O} \) is \( \mathbf{A} \)-invariant.

Solution: If \( \mathbf{x}\in\mathscr{N} \), then \( \mathbf{C}\mathbf{A}^k\mathbf{x}=\mathbf{0} \) for \( k=0,\ldots,n-1 \). For \( \mathbf{A}\mathbf{x} \), the observability rows give \( \mathbf{C}\mathbf{A}^k(\mathbf{A}\mathbf{x}) =\mathbf{C}\mathbf{A}^{k+1}\mathbf{x} \). The last condition is again zero by Cayley-Hamilton.

Problem 3: Let \( n=8 \), \( \dim\mathscr{R}=5 \), \( \dim\mathscr{N}=4 \), and \( \dim(\mathscr{R}\cap\mathscr{N})=2 \). Find the four block dimensions.

Solution:

\[ n_2=2,\qquad n_1=5-2=3,\qquad n_4=4-2=2,\qquad n_3=8-5-4+2=1. \]

Problem 4: Explain why the controllable-unobservable block may be input-reachable but absent from the transfer matrix.

Solution: The input may excite \( \mathbf{z}_2 \) through \( \mathbf{B}_2 \). However, \( \mathbf{C}_2=\mathbf{0} \), and \( \mathbf{A}_{12}=\mathbf{0} \), so this block cannot create output directly or indirectly through the controllable-observable block.

Problem 5: Compute the transfer function of the worked example.

Solution:

\[ G(s)=\mathbf{C}_1(s\mathbf{I}-\mathbf{A}_{11})^{-1}\mathbf{B}_1 =2(s+1)^{-1}=\frac{2}{s+1}. \]

Kalman decomposition reveals the internal anatomy of a realization: controllable-observable, controllable-unobservable, uncontrollable-observable, and uncontrollable-unobservable dynamics. The controllable-observable block is the minimal input-output subsystem. This result leads directly to Lesson 4, where nonminimal states are removed to construct a minimal realization.