Singular values of a matrix written in controllable canonical form Let the following equation represent a stable (marginally) dynamical system in discrete time domain
\begin{equation}
\mathbf{x}_{k+1} = \mathbf{A}\mathbf{x}_k + \mathbf{B}\mathbf{u}_k
\end{equation}
where $\mathbf{A}\in\mathbb{R}^{m\times m}$, $\mathbf{x}_k\in\mathbb{R}^m$, $\mathbf{B}\in\mathbb{R}^{m\times n}$ and $\mathbf{u}_k\in\mathbb{R}^{n}$. $\mathbf{A}$ is written in controllable canonical form. Since the system is stable (marginally), all the eigenvalues of $\mathbf{A}$ lie inside or on the unit circle. 

Given the system is stable, is it correct to say that at least $m-1$ singular values of $\mathbf{A}$ are less than or equal to unity and at most one singular value can be greater than $1$?

 A: A controller canonical form in control theory is just a companion matrix in linear algebra. It is known in linear algebra literature that for any $m\times m$ companion matrix (with $m\ge2$)
$$A=\pmatrix{
0 & 1 & 0 & \cdots & 0\\
0 & 0 & 1 & \cdots & 0\\
\vdots & \vdots & \vdots & \ddots & \vdots\\
0 & 0 & 0 & \cdots & 1\\
-a_0 & -a_1 & -a_2 & \cdots & -a_{m-1}},
$$
its middle $m-2$ singular values are given by $\sigma_2=\cdots=\sigma_{m-1}=1$ and its largest and smallest singular values are given by
$$
\sigma_1,\sigma_m=\sqrt{\frac{1+\|\mathbf a\|^2 \pm \sqrt{(1+\|\mathbf a\|^2)^2 - 4a_0^2\,}}2},
$$
where $\|\mathbf a\|$ is the Euclidean norm of $\mathbf a=(a_0,a_1,\ldots,a_{m-1})^T$. (See, e.g. Horn and Johnson, Matrix Analysis (2/e), p.197, 3.3.P13 or Alan Laub, Matrix Analysis for Scientists and Engineers, p.106, theorem 10.38.)
In your case, if the spectral radius of $A$ is at most $1$, we have $|a_0|\le1$. Therefore the only situation where $\sigma_1>1>\sigma_m$ fails to hold is that $|a_0|=1$ and $a_1=\cdots=a_{m-1}=0$.
A: Your claim is correct. In fact $\sigma=1$ is a singular value with multiplicity at least $m-2$.

If $\mathbf{A}$ has the controller canonical form
$$\mathbf{A}=\left[\matrix{0 & 1 & 0 & \cdots & 0\\ 0 & 0 & 1 & \cdots & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & 0 & \cdots & 1\\ -a_0 & -a_1 & -a_2 & \cdots & -a_{m-1}}\right]\qquad \qquad (1)$$ 
then the characteristic polynomial of $\mathbf{A}$  is
$$p_{\mathbf{A}}(\lambda)=\det(\lambda I-\mathbf{A})=\lambda^m+a_{m-1}\lambda^{m-1}+\cdots+a_0$$
It  holds true that  $a_0=(-1)^m\prod_{i=1}^m{\lambda_i(\mathbf{A})}$ and since $|\lambda_i(\mathbf{A})|<1$ for all $i=1,\cdots,m$ we directly deduce that $|a_0|<1$.
Then we can write
$$\mathbf{A^TA}=\left[\matrix{a_0^2 & a_0a_1 & a_0a_2 & \cdots & a_0a_{m-1}\\ a_0a_1 & 1+a_1^2 & a_1a_2 & \cdots & a_1a_{m-1}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ a_0a_{m-2} & a_1a_{m-2} & a_2a_{m-2} & \cdots & a_{m-2}a_{m-1}\\a_0a_{m-1} & a_1a_{m-1} & a_2a_{m-1} & \cdots & 1+a_{m-1}^2}\right]$$ 
which can be decomposed as

$$\mathbf{A^TA}=\mathbf{N}+\mathbf{a}\mathbf{a}^T$$
  with
  $$\mathbf{N}=diag\{0,1,1,\cdots,1\}\:, \quad\mathbf{a}=\left[\matrix{a_0 & a_1 & \cdots & a_{m-1}}\right]^T$$

If $\lambda$ is an eigenvalue of $\mathbf{A^TA}$ and $\mathbf{v}=\left[\matrix{v_1 & v_2 & \cdots & v_{m}}\right]^T$ the respective right eigenvector then
$$\lambda\mathbf{v}=\mathbf{A^TAv}=\mathbf{Nv}+\mathbf{a}\mathbf{a}^T\mathbf{v}$$
or equivalently

$$\lambda v_1  =(\mathbf{a}^T\mathbf{v})a_0\qquad \qquad \qquad \qquad \qquad \qquad\quad (2)\\  
 (\lambda-1) v_i=(\mathbf{a}^T\mathbf{v})a_{i-1}\:,\quad (i=2,\cdots,m) \qquad (3)$$

Define now the reduced vectors 
$$\mathbf{a_r}=\left[\matrix{a_1 & a_2 &\cdots & a_{m-1}}\right]^T\:,\quad \mathbf{v_r}=\left[\matrix{v_2 & v_3 &\cdots & v_{m}}\right]^T$$ and the vector $\mathbf{e_i}$ which denotes the $i$-th column of the identity matrix ($i=1,2,\cdots,m$).

Case 1: $\quad\mathbf{a_r}=0$
Then from (2), (3) we have 
$$\lambda v_1=a_0^2v_1\\ (\lambda-1)v_i=0\:,\quad (i=2,\cdots,m)$$
which mean that there are exactly 2 eigenvalues:

$\lambda_1=1$ with multiplicity $\mu_1=m-1$ and independent
  eigenvectors $\mathbf{e_2},\cdots, \mathbf{e_m}$ and a simple
  eigenvalue $\lambda_2=a_0^2<1$ with eigenvector $\mathbf{e_1}$.

Thus for Case 1 the maximum singular value of $\mathbf{A}$ is 1.

Case 2: $\quad\mathbf{a_r}\neq 0$
Then $\dim \ker(\mathbf{a_r}^T)=m-2$ (there are exactly $m-2$ independent vectors normal to $\mathbf{a_r}$) and from (2), (3) $\mathbf{A^TA}$ has an eigenvalue 

$\lambda_1=1$ with multiplicity $\mu_1=m-2$ and independent
  eigenvectors of the form $\left[\matrix{0 & \mathbf{v_r}^T}\right]^T$ with $\mathbf{v_r}\in \ker(\mathbf{a_r}^T)$.

We will examine now the behavior of the other 2 eigenvalues of $\mathbf{A^TA}$.

Case 2.A:$\quad a_0=0$
From (2), (3) we have

$$\lambda v_1  =0\qquad \qquad \qquad \qquad \\  
 (\lambda-1) \mathbf{v_r}=(\mathbf{a_r}^T\mathbf{v_r})\mathbf{a_{r}} \qquad $$

and therefore the other 2 eigenvalues are

$\lambda_2=0$ with eigenvector $\mathbf{e_1}$ and
  $\lambda_3=1+\mathbf{a_r^Ta_r}$ with eigenvector $\left[\matrix{0&
 \frac{\mathbf{a_r^T}}{\|\mathbf{a_r}\|}}\right]^T$

Thus, for Case 2.A there is exactly one singular value of $\mathbf{A}$ larger than one namely $\sqrt{1+\mathbf{a_r^Ta_r}}$.

Case 2.B:$\quad a_0\neq 0$
From (1) we can easily calculate that 
$\det(\mathbf{A})=(-1)^ma_0$ and therefore
$$\det(\mathbf{A^TA})=(\det(\mathbf{A}))^2=a_0^2$$
i.e. $0\neq \det(\mathbf{A^TA})<1$. Thus,

$$0\neq
 \det(\mathbf{A^TA})=\lambda_1^{m-2}\lambda_2\lambda_3=\lambda_2\lambda_3=a_0^2<1$$

and therefore at most one of the eigenvalues $\lambda_2,\lambda_3$ can be larger that one.
