Why is there this relationship between quaternions and Pauli matrices? I've just started studying quantum mechanics, and I've come across this correlation between Pauli matrices ($\sigma_i$) and quaternions which I can't grasp: namely, that $i\sigma_1$, $i\sigma_2$ and $i\sigma_3$, along with the 2x2 identity matrix $I$, correspond identically to the four 2x2 matrix representation of unit quaternions.
My first guess was that this should have something to do with quaternions being useful for representing orientations and rotations of objects in three dimensions and Pauli matrices being related to the three spatial components of spin, but I didn't really know how to put together those two ideas. Google wasn't much help either: the relation is mentioned, for instance, in this Wikipedia article, but no further explanation is given.
Although I suspect there is no direct answer to this question, I would appreciate if someone could enlighten me on the subject. In particular, what is the role of the $i$ factor?
 A: QMechanics has given you the straight answer in terms of group isomosphisms. So please go with that, but in case you go further into mathematics of quarterions and their applications in physics you will find many twists and turns on this subject.
I personally found the book On Quaternions and Octonions : Conway, Smith (2003) finally gave me some clarity on this whole subject. I will summarise a few key points. Apologies if this goes a bit further than your original question.
Qauternions are part of a series of division algebras used by mathematicians. They only occur in dimensions which are powers of 2 but only up to 8, namely:
1. Real numbers
2. Complex numbers
4. Quaterions
8. Octonians
As you should know, unit complex numbers are related to rotations in 2 dimensions - you might expect this to be part of a pattern and indeed it is roughly stated (from page 89):


*

*multiplication by unit complex numbers generate rotations in 2 dimensions 

*multiplication by unit quaternions generate rotations in 4
dimensions (not 3 dimensions - see below!) 

*multiplication by unit octonions generate rotations in 8 dimensions


The subtlety (which links back to your question) is that in 4 dimensions there are 6 rotations (you will know this if you have studied special relativity) so you actually need 2 copies of the quaternions. If you restrict to just one copy you get back to 3 dimension rotations.
In summary:
3d rotations: one copy of unit quaternions relate to pauli matrices
4d rotations: two copy of unit quaternions relate to 2 copies of pauli matrices
In group language:
Spin(3) = SU(2)           (3 dimensions)
Spin(4) = SU(2) x SU(2)   (6 dimensions)
As 3 and 4 dimensions are the two most important to physics this comes up in many guises in Quantum Physics. 
A: Summary: The Pauli matrices span the vector space of $2\times2$  traceless Hermitian matrices, the unit quaternions span the Lie algebra $2\times2$ traceless skew-Hermitian matrices (whence the $i$ factor), the latter being the Lie algebra of the Lie group of rotations and that group's universal cover.
As a particular basis for the lie algebra, the unit quaternions' defining relationships follow (1) from their being orthonormal with respect to the unique (up to a scaling constant) billinear form (i.e. here inner product) that is invariant under the action of $SO(3)$ (this gives the commutator relationships) and (2) requiring a $2\times2$ matrix representation of the Lie algebra. The Pauli matrices are then these elements divided by $i$, and you can think of the factor $i$ as being what is needed to make the Hermitian matrices close under a Lie bracket. 

The Pauli matrices span the vector space of traceless, $2\times 2$ Hermitian matrices and the quaternion units span the vector space of traceless, skew-Hermitian matrices, when we think of the faithful matrix representation of the vector space of so-called pure quaternions of the form $\alpha\,\mathbf{i}+\beta\,\mathbf{j}+\gamma\,\mathbf{k}$ where $\alpha,\,\beta,\,\gamma\in\mathbb{R}$. Therefore, the Pauli matrices must be linear superpositions of the scalar $i\in\mathbb{C}$ times the quaternion units $\mathbf{i},\, \mathbf{j},\,\mathbf{k}$ and indeed are related by a simple scaling by $i$ as you have found and as QMechanic's Answer elaborates with much detail.
The question is then: what has all this to do with rotation and spin?
To answer this, let's work in the other direction from the other answers by beginning with the group $SO(3)$ of rotations and seeing how the quaternions and Pauis follow from the geometry.
$SO(3)$, the group of proper, homogeneous isometries of Euclidean 3-space obviously captures the geometry of rotations. From this starting point, we need the following Lie group theoretic facts to go forward:


*

*$SU(2)$ is the double cover, indeed the universal cover of $SO(3)$;

*$SO(3)$ is the image of $SU(2)$ under the Adjoint representation $\mathrm{Ad}:SU(2)\to SO(3)$ whereby $SU(2)$ acts on its own Lie algebra $\mathfrak{su}(2)\cong\mathfrak{so}(3)$ by the spinor map. Thus $SU(2)$ acts on $\mathfrak{su}(2)$ by rigid rotations;

*$SO(3)$ and $SU(2)$ are simple Lie groups, containing no normal proper Lie subgroups of dimension 1 or greater and thus $\mathfrak{su}(2)\cong\mathfrak{so}(3)$ is a simple Lie algebra;

*The Killing form for the Lie algebra $\mathfrak{g}$ of a simple Lie group $\mathfrak{G}$ is the unique billinear form, to within a real scaling constant, that is $\mathrm{Ad}$-invariant. That is, if $X,\,Y \in\mathfrak{g}$ and $F(X,\,Y)$ is the value of a billinear form defined on the Lie algebra, and if further $F(X,\,Y)=F(\mathrm{Ad}(\gamma)\,X, \, \mathrm{Ad}(\gamma)\,Y);\forall X,\,Y\in\mathfrak{g},\,\forall \gamma\in\mathfrak{G}$, then $F$ is the Killing form, to within a real scaling constant.
So now the question is, what's the simplest basis for the Lie algebra $\mathfrak{su}(2)$ that exponentiates to our rotation group and its double cover?
$SU(2)$ acts, through the Adjoint representation, on $\mathfrak{su}(2)$ by rigid rotations, so let's define an inner product on $\mathfrak{su}(2)$ that is left invariant by the group action on the algebra. This inner product will then yield the correct angles between vectors in $\mathfrak{su}(2)$ since this is what is left invariant by rotations.
But we know $SU(2),\,SO(3)$ are simple, so there is only one billinear form that is invariant in this way to within a scaling constant: the Killing form.
So now, it is obviously convenient to construct a basis that is orthonormal with respect to the Killing form.
With a little bit of work, one can show that such an orthonormal basis for $\mathfrak{su}(2)$  must fulfill $[\mathbf{i},\,\mathbf{j}]=2\,\mathbf{k}$, $[\mathbf{k},\,\mathbf{i}]=2\,\mathbf{j}$, $[\mathbf{j},\,\mathbf{k}]=2\,\mathbf{i}$. Furthermore, if we use $2\times2$ unitary matrices to represent $SU(2)$, then such an orthonormal basis for $\mathfrak{su}(2)$ must also fulfill the famous "Broughham bridge" relationships (equation (1) in QMechanic's Answer).
So now our quaternion units are defined, to within a similarity transformation. One obvious choice is the one that gives the relationship between the Paulis and the quaternion units that you cite.
A: *

*At the level of formulas, the three quaternionic units $i_a$,  $a\in~\{1,2,3\}$, in $\mathbb{H}\cong \mathbb{R}^4$ satisfy
$$i_a i_b ~=~ -\delta_{ab}  + \sum_{c=1}^3\varepsilon_{abc} i_c, \qquad\qquad a,b~\in~\{1,2,3\}, \tag{1}$$
while the three Pauli matrices $\sigma_a \in {\rm Mat}_{2\times 2}(\mathbb{C})$,  $a\in~\{1,2,3\}$, $\mathbb{C}=\mathbb{R}+\mathrm{i}\mathbb{R}$, satisfy
$$\sigma_a \sigma_b ~=~ \delta_{ab} {\bf 1}_{2\times 2} + \mathrm{i}\sum_{c=1}^3\varepsilon_{abc} \sigma_c\quad\Leftrightarrow \quad \sigma_{4-a} \sigma_{4-b} ~=~ \delta_{ab} {\bf 1}_{2\times 2} - \mathrm{i}\sum_{c=1}^3\varepsilon_{abc} \sigma_{4-c}, $$
$$ \qquad\qquad a,b~\in~\{1,2,3\},\tag{2}$$
with complex unit $\mathrm{i}\in\mathbb{C}.$
In other words, we evidently have an $\mathbb{R}$-algebra monomorphism
$$\Phi:~~\mathbb{H}~~\longrightarrow ~~{\rm Mat}_{2\times 2}(\mathbb{C}).\tag{3}$$
by extending the definition 
$$\Phi(1)~=~{\bf 1}_{2\times 2},\qquad \Phi(i_a)~=~\mathrm{i}\sigma_{4-a}, \qquad\qquad a~\in~\{1,2,3\},\tag{4}$$
via $\mathbb{R}$-linearity. This observation essentially answers OP title question (v2). 

*However OP's question touches upon many beautiful and useful mathematical facts about Lie groups and Lie algebras, some of which we would like to mention. The image of the $\mathbb{R}$-algebra monomorphism (3) is
$$\Phi(\mathbb{H})
~=~ \left\{\left. \begin{pmatrix} \alpha & \beta \cr -\bar{\beta} & \bar{\alpha} \end{pmatrix}\in {\rm Mat}_{2\times 2}(\mathbb{C}) \right| \alpha,\beta \in\mathbb{C}\right\}$$
$$~=~ \left\{ M\in {\rm Mat}_{2\times 2}(\mathbb{C}) \left| \overline{M} \sigma_2=\sigma_2 M\right. \right\}.\tag{5}$$
Let us for the rest of this answer identify $\mathrm{i}=i_1$. Then the $\mathbb{R}$-algebra monomorphism (3) becomes 
$$ \mathbb{C}+\mathbb{C}i_2~=~\mathbb{H}~\ni~x=x^0+\sum_{a=1}^3 i_a x^a ~=~\alpha+\beta i_2$$
$$~~\stackrel{\Phi}{\mapsto}~~ \begin{pmatrix} \alpha & \beta \cr -\bar{\beta} & \bar{\alpha} \end{pmatrix} ~=~ x^0{\bf 1}_{2\times 2}+\mathrm{i}\sum_{a=1}^3 x^a \sigma_{4-a}~\in~ {\rm Mat}_{2\times 2}(\mathbb{C}),$$
$$ \alpha~=~x^0+\mathrm{i}x^1~\in~\mathbb{C},\qquad  \beta~=~x^2+\mathrm{i}x^3~\in~\mathbb{C},\qquad x^0, x^1, x^2, x^3~\in~\mathbb{R}.\tag{6}$$

*One may show that $\Phi$ is a star algebra monomorphism, i.e. the Hermitian conjugated matrix satisfies
$$ \Phi(x)^{\dagger}~=~\Phi(\bar{x}), \qquad x~\in~\mathbb{H}. \tag{7}$$
Moreover, the determinant becomes the quaternionic norm square
$$\det \Phi(x)~=~ |\alpha|^2+|\beta|^2~=~\sum_{\mu=0}^3 (x^{\mu})^2 ~=~|x|^2, \qquad x~\in~\mathbb{H}.\tag{8}$$
Let us for completeness mention that the transposed matrix satisfies 
$$\Phi(x)^t~=~\Phi(x|_{x^2\to-x^2})~=~ \Phi(-j\bar{x}j), \qquad x~\in~\mathbb{H}. \tag{9} $$ 

*Consider the Lie group of quaternionic units, which is also the Lie group 
$$U(1,\mathbb{H})~:=~\{x\in\mathbb{H}\mid |x|=1 \} \tag{10}$$ 
of unitary $1\times 1$ matrices with quaternionic entries. Eqs. (7) and (8) imply that the restriction $$\Phi_|:~U(1,\mathbb{H})~~\stackrel{\cong}{\longrightarrow}~~ SU(2)~:=~\{g\in {\rm Mat}_{2\times 2}(\mathbb{C})\mid g^{\dagger}g={\bf 1}_{2\times 2},~\det g = 1  \} $$
$$~=~\left\{\left. \begin{pmatrix} \alpha & \beta \cr -\bar{\beta} & \bar{\alpha} \end{pmatrix} \in {\rm Mat}_{2\times 2}(\mathbb{C}) \right| \alpha, \beta\in\mathbb{C}, |\alpha|^2+|\beta|^2=1\right\}\tag{11}$$ 
of the monomorphism (3) is a Lie group isomorphism. In other words, we have shown that

$$ U(1,\mathbb{H})~\cong~SU(2).\tag{12}$$


*Consider the corresponding Lie algebra of imaginary quaternionic number
$$ {\rm Im}\mathbb{H}~:=~\{x\in\mathbb{H}\mid x^0=0 \}~\cong~\mathbb{R}^3 \tag{13}$$
endowed with the commutator Lie bracket. [This is (twice) the usual 3D vector cross product in disguise.] The corresponding Lie algebra isomorphism is
$$\begin{align}\Phi_|:~{\rm Im}\mathbb{H}~~\stackrel{\cong}{\longrightarrow}~~ 
su(2)~:=~&\{m\in {\rm Mat}_{2\times 2}(\mathbb{C})\mid m^{\dagger}=-m \}\cr
~=~&\mathrm{i}~{\rm span}_{\mathbb{R}}(\sigma_1,\sigma_2,\sigma_3),\end{align}\tag{14}$$ 
which brings us back to the Pauli matrices. In other words, we have shown that

$$ {\rm Im}\mathbb{H}~\cong~su(2).\tag{15}$$


*It is now also easy to make contact to the left and right Weyl spinor representations  in 4D spacetime $\mathbb{H}\cong \mathbb{R}^4$ endowed with the quaternionic norm $|\cdot|$, which has positive definite Euclidean (as opposed to Minkowski) signature, although we shall only be sketchy here. See also e.g. this Phys.SE post. 

Firstly, $U(1,\mathbb{H})\times U(1,\mathbb{H})$  is (the double cover of) the special orthogonal group $SO(4,\mathbb{R})$. 

The group representation
$$\rho: U(1,\mathbb{H}) \times U(1,\mathbb{H}) \quad\to\quad SO(\mathbb{H},\mathbb{R})~\cong~ SO(4,\mathbb{R}) \tag{16}$$
is given by
$$\rho(q_L,q_R)x~=~q_Lx\bar{q}_R, \qquad q_L,q_R~\in~U(1,\mathbb{H}), \qquad x~\in~\mathbb{H}. \tag{17}$$ 
The crucial point is that the group action (17) preserves the norm, and hence represents orthogonal transformations. See also this math.SE question.

Secondly, $U(1,\mathbb{H})\cong SU(2)$  is (the double cover of) the special orthogonal group $SO({\rm Im}\mathbb{H},\mathbb{R})\cong SO(3,\mathbb{R})$. 

This follows via a diagonal restriction $q_L=q_R$ in eq. (17).
A: The $σ$ matrices generate a Clifford algebra, like the Dirac $γ$ matrices. In this case it is the three-dimensional algebra of signature $+{+}+$.
In odd dimensions, the volume element/pseudoscalar of a Clifford algebra (in this case $\mathbf{xyz}$ or $\mathbf e_1 \mathbf e_2 \mathbf e_3$ or $σ_xσ_yσ_z$) commutes with everything. If in addition it squares to $-1$, which in this signature it does, then it can be identified with the complex $i$. Note that $σ_xσ_yσ_z = \mathrm{diag}(i,i)$. The $σ$ matrices therefore generate the real Clifford algebra only, and $i$ represents the Clifford pseudoscalar and not the separate $i$ that would appear in a complexified algebra. This is different from the Dirac case: in 3+1 dimensions, while the pseudoscalar still squares to $-1$, it anticommutes with the odd part of the algebra, and so can't be identified with $i$.
The even subalgebra of any Clifford algebra contains a representation of the spin group (and the full algebra contains a representation of the pin group). In signature $+{+}+$ and $-{-}-$, the even subalgebra is isomorphic to the quaternions, and elements of norm $1$ represent elements of $\mathrm{Spin}(3)$ in the same way that unit quaternions do.
The basis quaternions can be written as pairwise products of $σ$ matrices or as products of single $σ$ matrices with the pseudoscalar. But note that you can't take $\mathbf i, \mathbf j, \mathbf k = iσ_x, iσ_y, iσ_z$, because then you get $\mathbf{ijk}=+1$, which is fine in the abstract but isn't what Hamilton chose. You need $\mathbf i, \mathbf j, \mathbf k = -iσ_x, -iσ_y, -iσ_z$ instead. (You could also take any odd permutation of $iσ_x, iσ_y, iσ_z$, but that makes no geometric sense if you're associating the $σ$ matrices and basis quaternions with coordinate axes, as you should.)

You can alternately identify the $σ$ matrices with the generating set $\pm\mathbf{tx}, \pm\mathbf{ty}, \pm\mathbf{tz}$ of the even, real Clifford algebra of 3+1 dimensions, where the sign should be chosen so that $σ_xσ_yσ_z$ is the pseudoscalar. Because this is only the even subalgebra, the pseudoscalar is central as before and it can be (and is) identified with $i$. The quaternions appear as the subalgebra of rotations in the hypersurface normal to $\mathbf t$. The 6-dimensional subalgebra isomorphic to $\mathrm{Spin}(3,1)$ turns out to be $\mathrm{SL}(2,\mathbb C)\subset \mathrm M_2(\mathbb C)$.
A: There is a simple answer to your question. First, in $\mathrm{Cl3}$ (geometric algebra of 3D Euclidean vector space), the elements of the even part of the algebra are just quaternions.
Second, the simplest representation of orthonormal vectors in $\mathrm{Cl3}$ are the Pauli matrices. In addition, we have a geometric meaning of both the Pauli matrices and quaternions.
A: All the answers above are very good and illuminating mathematically. I like to add a simple physical picture. See this article and quote below; http://wwwf.imperial.ac.uk/~jdg/nonlinquat.pdf 
‘’A quaternion can be constructed from a scalar s and a 3-vector r by forming the tetrad q = [s,r] that is defined by q = [s, r] = sI − r • σ,  where r• σ = Sum(ri σi); i=1-3, and I is the 2 × 2 unit matrix. {σ1, σ2, σ3} are the Pauli spin matrices σ1 = [0 1; i 0] , σ2 = [0 1; −1 0] , σ3 =[ i 0; 0 –i] , that obey the relations;  σiσj=−δij I – εijk σk . A multiplication rule between two tetrads q1 = [s1, r1] and q2 = [s2, r2] can easily be determined from these properties q1 * q2 = [s1s2 − r1 • r2, s1r2 + s2r1 + r1 × r2].’’
The article is in fluids, so what is the physics connecting the two.
Radiation or fields in general follow hyperbolic type equations- the wave equation. This equation represents advancing in space and time and doesn't allow any back flow/vortices/curvature/spin. To be able to have flow in the reverse direction we need an elliptic type equation. This can represent vortices and spin. Dirac equation is elliptic while Klein-Gordon is hyperbolic.
This is where we need the rotation matrices/the Pauli matrices or quaternions. 
In the flow of fluids, the wall does the job of deflecting the flow and creating vortices and turbulence. In the present case, it is done by the imaginary 'i', since multiplying by i rotates a quantity by 90 degrees as we know.   
A: For a very easy and quick explanation, use a search engine to find the paper titled "Imaginary Numbers are not Real". Eight pages in you get your answer.
Pauli matrices are mathematically and physically obtuse, a poor way to teach anything, let alone quantum mechanics. In a nutshell, quaternions are the proper way to think about spin in 3D. Imagine two objects spinning in two separate planes perpendicular to each other in that higher dimensional space. Let one "multiply" the other. Literally the first element in the operation, acts on the second, moving it into the 3rd plane that defines the space while maintaining its spin orientation. The operation is anti-commutative: reverse the order of the elements being multiplied, and you end up in that 3rd plane, but with spin in the opposite direction. 
That is how spinning particles interact. It's a function of the geometry of space, not a property of the particles. 
