This is in fact an interesting question. It involves some deep aspect about physics, symmetry, Lie groups and Lie algebra. There are many times when you can actually derive the Hermitian operators as physical observables without eigenvalues or even eigenvectors. You just need an appropriate coordinate system, some physical quantity, and symmetry.
Starting Point
When we have some physical quantity represented by coordinates, one natural question is: how to express the same physical quantity when we change our coordinates?
Translational Operator
For example, if we have a one-dimensional wavefunction $\psi(x)$, how would we express it if we move our origin to the left by $a$? We just have the new wavefunction (the wavefunciton in the new coordinate) $\tilde{\psi}(x)=\psi(x-a)$. The answer is simple. However, if we go further from this, expressing
$$\tilde{\psi}(x) = \psi(x-a) = \sum_{k=0}^{+\infty}{{1 \over k!}(-a)^k{\partial^k \over \partial x^k}\psi(x)} = e^{-a{\partial \over \partial x}}\psi(x)$$
We have the translational operator $\hat{T}(a)=e^{-a{\partial \over \partial x}}=e^{-ia{1 \over i}{\partial \over \partial x}}$, which is a Lie group parametrized by $a \in \mathbb{R}$. The generator of the Lie group is $-{\partial \over \partial x}$, which is the basis of the Lie algebra. You probably recognize ${1 \over i}{\partial \over \partial x}={\hat{p}_x \over \hbar}$. This is not a special case, and you can generalize the result to the operator $\hat{p}\cdot{\hat{n}}$ where $-\hat{n}$ is the direction we move the origin. Also, if you have some time-independent Hamiltonian $\hat{H}$, and translate the time coordinate to the left by $t$, you will have the Lie group $e^{i t {\hat{H} \over \hbar}}$ which connects $\tilde{\psi}(t')=\psi(t'-t)=e^{i t {\hat{H} \over \hbar}}\psi(t')$, and the basis of the Lie algebra $i{\hat{H} \over \hbar}$.
Spin-${1 \over 2}$
What about spin? Recall that spin is the internal angular momentum which an elementary particle has. So it is reasonable to see how spin states change if we rotate, instead of translating, the coordinate system. The spin-${1 \over 2}$ system is the most simple one to start with, also the one included in your question. Let us put down the axiom about spin-${1 \over 2}$ here:
- If we set up a coordinate system, there is a spin-up state $|\uparrow\rangle$ and a spin-down state $|\downarrow\rangle$. If we measure the angular momenta of these states, the only non-zero components are in the $z$-axis with the value of $|\uparrow\rangle$ higher than $|\downarrow\rangle$.
- $|\uparrow\rangle$ and $|\downarrow\rangle$ compose a orthonormal basis: Any spin state $|\psi\rangle$ is a superposition of spin-up and spin-down, so we have $|\psi\rangle=a|\uparrow\rangle+b|\downarrow\rangle=\begin{pmatrix}{a \\ b}\end{pmatrix}$ where $a,b \in \mathbb{C}$ and $|a|^2+|b|^2=1$. Also, the converse is true. If $a,b \in \mathbb{C}$ and $|a|^2+|b|^2=1$, we have the spin state $|\psi\rangle=a|\uparrow\rangle+b|\downarrow\rangle=\begin{pmatrix}{a \\ b}\end{pmatrix}$.
Let us have some spin state $|\psi\rangle$ and rotate the coordinate system. As results, we have the new spin state (the spin state expressed in the new coordinate) $|\psi'\rangle=\begin{pmatrix}{a' \\ b'}\end{pmatrix}=a'|\uparrow'\rangle+b'|\downarrow'\rangle$ where $|a'|^2+|b'|^2=1$ and $|\uparrow'\rangle$ and $|\downarrow'\rangle$ are the spin-up and spin-down state in the new coordinate system. This gives us a $2\times2$ rotation matrix $U$ which satisfies $\begin{pmatrix}{a' \\ b'}\end{pmatrix}=U\begin{pmatrix}{a \\ b}\end{pmatrix}$. Since $a'$ and $b'$ are arbitrary, we have a set of $U$, which is
$$\text{SU}(2)=\bigg\{\begin{pmatrix}u_1 & u_2 \\ -u^{\dagger}_2 & u^{\dagger}_1\end{pmatrix}\bigg|u_1,u_2 \in \mathbb{C}, |u_1|^2+|u_2|^2=1\bigg\}$$
So the problem becomes: How are we going to relate the rotation of our three dimensional space, which belongs to $\text{SO}(3)$, to the element of $\text{SU}(2)$? Luckily, we have the mathematical statement that $\text{SO}(3) \cong \text{SU}(2) / \{\pm I_2\}$ where $I_2$ is the $2\times2$ identiy matrix. This means if we rotate our coordinate system around some axis $\hat{n}=(n_x,n_y,n_z)$ by $-\alpha$ degrees, the corresponding element $U$ in $SU(2)$ has to be
\begin{align}
e^{-i\alpha{1 \over 2}\vec{\sigma} \cdot \hat{n}} & = e^{-i\alpha{1 \over 2}(\sigma_xn_x+\sigma_yn_y+\sigma_zn_z)} \\
& = e^{-i\alpha{1 \over 2}\bigg(n_x\begin{pmatrix}0 & 1 \\ 1 & 0\end{pmatrix}+n_y\begin{pmatrix}0 & -i \\ i & 0\end{pmatrix}+n_z\begin{pmatrix}1 & 0 \\ 0 & -1\end{pmatrix}\bigg)}
\end{align}
The generators of the Lie group $\text{SU}(2)$ are $-i{\hat{S}_x \over \hbar}$, $-i{\hat{S}_y \over \hbar}$, and $-i{\hat{S}_z \over \hbar}$, which are the basis of the Lie algebra $\text{su}(2)$.
At this point, you may already summarize a way to derive Hermitian operators:
- We unitarily transform the coordinate system given some physical quantity,
- and express the transformations in terms of unitary Lie groups. The generators are necessarily skew-Hermitian but give us physical Hermitian operators after multiplication of $i$.
However, why does this work?
Symmetry
Let us go back to our translational operators $\hat{T}(a)$, and assume we do not know the form of momentum operator $\hat{p}_x$. However, there is a property we require $\hat{p}_x$ to have:
$$\int{\psi^{\dagger}(x-a)\hat{p}_x\psi(x-a)dx}=\text{const}$$
for any wavefuntion $\psi$ and $a$. The requirement is reasonable since we do not think the momentum of a wavefunction varies if we just translate our coordinate. This means the quantity $\int{\psi^{\dagger}(x-a)\hat{p}_x\psi(x-a)dx}$ has the translational symmetry, and it implies
\begin{align}
{d \over da}\int{\psi^{\dagger}(x-a)\hat{p}_x\psi(x-a)dx} & = {d \over da}\int{(\hat{T}(a)\psi(x))^{\dagger}\hat{p}_x\hat{T}(a)\psi(x)dx} \\
& = {d \over da}\int{\psi(x)\big(\hat{T}(-a)\hat{p}_x\hat{T}(a)\big)\psi(x)dx} \\
& = \int{\psi(x){d \over da}\big(e^{a{\partial \over \partial x}}\hat{p}_xe^{-a{\partial \over \partial x}}\big)\psi(x)dx} = 0
\end{align}
Therefore, ${d \over da}\big(e^{a{\partial \over \partial x}}\hat{p}_xe^{-a{\partial \over \partial x}}\big)=0$, and $e^{a{\partial \over \partial x}}\hat{p}_xe^{-a{\partial \over \partial x}}=\hat{p}_x$. Through limiting $a \rightarrow 0$, we have
$$\hat{p}_x+a({\partial \over \partial x}\hat{p}_x-\hat{p}_x{\partial \over \partial x})+O(a^2)=\hat{p}_x$$
So this requires $\hat{p}_x$ to commute with ${\partial \over \partial x}$. This means we can restrict $\hat{p}_x$ to have the same eigenvectors as ${\partial \over \partial x}$. At least, it makes a good guess that $\hat{p}_x=ik{\partial \over \partial x}$ where $k \in \mathbb{R}$ (there are other reasons showing $\hat{p}_x={\hbar \over i}{\partial \over \partial x}$ for plane waves, which helps tying up the loose end).
The same argument can be applied to the spin matrices. Take $\hat{S}_z$ as an example. A symmetry which $\hat{S}_z$ must have is the rotational symmetry about $z$-axis since any rotation about the $z$-axis should not change the expectation value of the angular momentum along the $z$-axis for any spin state $|\psi\rangle$. Therefore, $\hat{S}_z$ must commute with $\begin{pmatrix}1 & 0 \\ 0 & -1\end{pmatrix}$, and if we further impose that $\hat{S}_z$ is traceless (this is reasonable since if we rotate the cooridnate about $x$-axis by $\pi$, the spin-up state $|\uparrow\rangle$ is parallel to $|\downarrow'\rangle$ in the new coordinate), we have $\hat{S}_z=k\begin{pmatrix}1 & 0 \\ 0 & -1\end{pmatrix}$.
