Canonical transformation for the "nearly diagonal" Hamiltonians The Hamiltonians in question are the Hamiltonians that have all the non-diagonal equal to zero except for the first row and the first column (assuming for simplicity that $U_i\in\mathbb{R}$ are real):
$$
H_2 = 
\begin{bmatrix}
\epsilon_0 & U_1 \\
U_1 & \epsilon_1
\end{bmatrix}
$$
$$
H_3 = 
\begin{bmatrix}
\epsilon_0 & U_1 & U_2\\
U_1 & \epsilon_1 & 0 \\
U_2 & 0 & \epsilon_2
\end{bmatrix},
$$
$$
H_4 = 
\begin{bmatrix}
\epsilon_0 & U_1 & U_2 & U_3\\
U_1 & \epsilon_1 & 0 & 0 \\
U_2 & 0 & \epsilon_2 & 0 \\
U_3 & 0 & 0 & \epsilon_3
\end{bmatrix},
$$
and so on.
Relevance to physics
Such Hamiltonians are quite ubiquitous, and many people could come with the examples typical for their field. I will mention just a few:

*

*Coupled sublattice sites (order $n$ of $H_n$ is typically 3 or 4, see also this question)

*A level coupled to a band (order $n$ is very large)

*Anderson impurity (with some complications added)

Solvability
Many particular cases of this Hamiltonian are solvable: e.g., when all the couplings are the same ($U_1=U_2=...=U_{n_1}$) and the number of levels is infinite (broad-band limit for a resonant level) or if $\epsilon_1=\epsilon_2=...=\epsilon_{n-1}\neq \epsilon_0$, when the characteristic equation can be reduced to a quadratic equation.
Moreover, in cases $n=3$, $n=4$ the characteristic equation is exactly solvable in principle, although the results may not look pretty.
Question
What is the general form of a canonical transformation for diagonalizing such Hamiltonians?
For $n=3$ and $n=4$ this would amount to a simply convenient parametrization, whereas for arbitrary $n$ it is a transformation that takes into account the greater simplicity of this Hamiltonians in comparison to a general $n\times n$ Hamiltonian.
To suggest a specific line of reasoning: the canonical transformation matrix for a real Hamiltonian is an orthogonal matrix with $n(n-1)/2$ independent parameters (see here, e.g.). However, the zeros in the Hamiltonian impose additional $(n-1)(n-2)/2$ constraints, thus leaving only $n-1$ independent parameters.
Update
I am looking for a closed analytical expression parametrized in terms of $n-1$ parameters, which could be unknown. For example, this could be $n-1$ eigenvalues - since we are not sure that they can be found analytically for $n>4$ (a definitive statement to this end would be also handy). I might be already providing an answer here...
 A: Disclaimer: Due to the disagreement (in the comments) about what should be considered an answer to my question I post this here: it refines the existing answer and raises more precise questions about what needs to be understood for a complete treatment of this problem.
Standard approach to eigenvalue problems includes

*

*Solving the characteristic equation for the eigenvalues, $$|H-\lambda I|=0,$$ where $I$ is the identity matrix.

*Solving the eigenvectors equations, given known eigenvalues,
$$
H\mathbf{x}=\lambda \mathbf{x}.
$$
The characteristic equation
Evaluating determinants is is trivial for this type of the Hamiltonians. In particular we obtain:
$n=2$
$$
|H-\lambda I| = \left|\begin{matrix} \epsilon_0-\lambda & U_1 \\ U_1 & \epsilon_1 - \lambda \end{matrix}\right|
=(\epsilon_0 - \lambda)(\epsilon_1-\lambda) - U_1^2 = 0
$$
$n=3$
$$
|H-\lambda I| = 
(\epsilon_0 - \lambda)(\epsilon_1-\lambda)(\epsilon_2-\lambda) - U_1^2(\epsilon_2-\lambda) - U_2^2(\epsilon_1-\lambda) = 0
$$
$n=3$
$$
|H-\lambda I| = 
(\epsilon_0 - \lambda)(\epsilon_1-\lambda)(\epsilon_2-\lambda)(\epsilon_3-\lambda) - U_1^2(\epsilon_2-\lambda)(\epsilon_3-\lambda) - U_2^2(\epsilon_1-\lambda)(\epsilon_3-\lambda) -U_3^2(\epsilon_1-\lambda)(\epsilon_2-\lambda) = 0
$$
One can then show by induction that for an arbitrary order we get:
$$
|H-\lambda I| = (\epsilon_0 - \lambda)\prod_{i=1}^{n-1} (\epsilon_i - \lambda) -  \sum_{i=1}^{n-1}U_i^2 \prod_{j=1, j\neq i}^{n-1} (\epsilon_i - \lambda) = \\
\left(\epsilon_0 - \lambda - \sum_{i=1}^{n-1}\frac{U_i^2}{\epsilon_i-\lambda}\right)\prod_{i=1}^{n-1} (\epsilon_i - \lambda)=0
$$
The last form of the determinant obviously makes sense only if the roots are not equal to any of $\epsilon_i, i>0$. It is also easily seen that, if any two energies are equal $\epsilon_i=\epsilon_j$, one can factor out $\epsilon_i-\lambda$, obtaining one root of the equation.
Q1: The non-trivial question on the level of the characteristic equation is: can it be solved or can it be shown that it cannot be solved for $n>0$. Note that while the latter statement is known to be true for a general polynomial equation, it is not clear a priori that this is also the case for this particular form.
Solving for eigenvectors
Assuming that the eigenvalues $\lambda_i, i=1,..,n$ are known, we can solve for the eigenvectors (Note that this is again a standard procedure for an arbitrary matrix Hamiltonian). For simplicity we assume that none of the roots is equal to any of the energies $\epsilon_j$, although we will show that the solution can be represented in a general form, suitable for a degenerate case.
The solution for this types of the Hamiltonians is trivial (formally one could call it a trivial case of Gaussian elimination). The eigenvector corresponding to eigenvalue $\lambda_i$ takes form:
$$
\mathbf{x}_i = A_i\begin{bmatrix} 1\\\frac{U_1}{\lambda_i - \epsilon_1}\\\frac{U_2}{\lambda_i - \epsilon_1}\\...\\\frac{U_{n-1}}{\lambda_i - \epsilon_{n-1}}\end{bmatrix},
$$
where
$$
A_i = \left[1 + \sum_{j=1}^{n-1}\frac{U_j^2}{(\lambda_i - \epsilon_j)^2}\right]^{-\frac{1}{2}}.
$$
Note that by replacing $1$ in the first element in the eigenvector with $\prod_{j=1}^{n-1}(\lambda_i-\epsilon_j)$ we could remove singularities from thz denominators, obtaining the answers that are also valid for a non-degenerate case.
The canonical transformation diagonalizing the Hamiltonian ($S^\dagger H S = \Lambda$, where $\Lambda_{ij} = \delta_{i,j}\lambda_i$) can be written as
$$
\mathbf{S} = 
\begin{bmatrix} 
A_1 & A_2 & ... & A_3\\
\frac{A_1 U_1}{\lambda_1 - \epsilon_1} & \frac{A_2 U_1}{\lambda_2 - \epsilon_1} & ... & \frac{A_{n} U_1}{\lambda_{n} - \epsilon_1}\\
\frac{A_1 U_2}{\lambda_1 - \epsilon_2} & \frac{A_2 U_2}{\lambda_2 - \epsilon_2} & ... & \frac{A_{n} U_2}{\lambda_{n} - \epsilon_2}\\
... & ... & ... & ... \\
\frac{A_1 U_{n-1}}{\lambda_1 - \epsilon_{n-1}} & \frac{A_2 U_{n-1}}{\lambda_2 - \epsilon_{n-1}} & ... & \frac{A_{n} U_{n-1}}{\lambda_{n} - \epsilon_{n-1}}
\end{bmatrix}.
$$
This canonical transformation is parametrized by $n$ eigenvalues $\lambda_i$, which can be considered as parameters: even if we fail to solve the characteristic equations the canonical transformation will have this form.
Q2: Contrary to what was suggested in the OP, the transformation diagonalizing an arbitrary n-by-n Hamiltonian is parametrizable in terms of its $n$ eigenvalues. So the question is about the special constrains and the number of independent parameters that we have in case of the simpler Hamiltonians considered in this question. Given the presence of the additional constrains - shouldn't there exist a relation between the eigenvalues, meaning that some cases with $n>4$ are actually exactly solvable?
Simpler form of the canonical transformation
The canonical transformation above can be written as
$$
\mathbf{S} = 
\begin{bmatrix}
1 & 0 & ... & 0 \\
0 & \frac{U_1}{U} & ... & 0 \\
... & ... & ... & ... \\
0 & 0 & ... & \frac{U_{n-1}}{U}
\end{bmatrix}
\begin{bmatrix} 
A & A & ... & A\\
\frac{A U}{\lambda_1 - \epsilon_1} & \frac{A U}{\lambda_2 - \epsilon_1} & ... & \frac{A U}{\lambda_{n} - \epsilon_1}\\
\frac{A U}{\lambda_1 - \epsilon_2} & \frac{A U}{\lambda_2 - \epsilon_2} & ... & \frac{A U}{\lambda_{n} - \epsilon_2}\\
... & ... & ... & ... \\
\frac{A U}{\lambda_1 - \epsilon_{n-1}} & \frac{A U}{\lambda_2 - \epsilon_{n-1}} & ... & \frac{A U}{\lambda_{n} - \epsilon_{n-1}}
\end{bmatrix}
\begin{bmatrix}
\frac{A_1}{A} & 0 & ... & 0 \\
0 & \frac{A_2}{A} & ... & 0 \\
... & ... & ... & ... \\
0 & 0 & ... & \frac{A_n}{A}
\end{bmatrix}
$$
This decomposition is more obvious, if we take $U=A=1$, but defining these factors as
$$
U=\left(\prod_{i=1}^{n-1}U_i\right)^{\frac{1}{n-1}},\\
A=\left(\prod_{i=1}^{n-1}A_i\right)^{\frac{1}{n}}
$$
preserves the dimensionality and assures that we are dealing with a product of three canonical transformations (three unitary/orthogonal matrices).
Q3: The ultimate holy grail here would be to have a representation of the canonical transformation that is manifestly unitary and "symmetric" in the way the well-known rotation matrices are for $n=2$ and $n=3$. At present these properties are hidden in the eigenvalues $\lambda_i$.
Other possible approaches

*

*Canonical transformation with an ad-hoc parametrization of the transformation matrix (e.g., using Cayley transform)

*Formal perturabtion expansion to all orders of magnitude, e.g., using the resolvent operator

*Using an auxhiliary transformation, as suggested by @Eli

