I have read about singular value decomposition (SVD). In almost all textbooks it is mentioned that it factorizes the matrix into three matrices with given specification. 

But what is the intuition behind splitting the matrix in such form? PCA and other algorithms for dimensionality reduction are intuitive in the sense that algorithm has nice visualization property but with SVD it is not the case. 

Write the SVD of matrix $X$ (real, $n\times p$) as
$$
   X = U D V^T
$$
where $U$ is $n\times p$, $D$ is diagonal $p\times p$ and $V^T$ is $p\times p$. In terms of the columns of the matrices $U$ and $V$ we can write
$X=\sum_{i=1}^p d_i u_i v_i^T$. That shows $X$ written as a sum of $p$ rank-1 matrices. What does a rank-1 matrix look like? Let's see:
$$
\begin{pmatrix} 1 \\ 2 \\ 3 \end{pmatrix} \begin{pmatrix} 4 & 5 & 6 \end{pmatrix} = \begin{pmatrix} 4 & 5 & 6 \\ 8 & 10 & 12 \\ 12 & 15 & 18 \end{pmatrix}
$$ The rows are proportional, and the columns are proportional.
Think now about $X$ as containing the grayscale values of a black-and-white image, each entry in the matrix representing one pixel.  For instance the following picture of a baboon:

Then read this image into R and get the matrix part of the resulting structure, maybe using the library pixmap.

If you want a step-by-step guide as to how to reproduce the results, you can find the code here.

Calculate the SVD:
    baboon.svd  <-  svd(bab) # May take some time

How can we think about this?  We get the $512 \times 512$ baboon image represented as a sum of $512$ simple images, with each one only showing vertical and horizontal structure, i.e. it is an image of vertical and horizontal stripes! So, the SVD of the baboon represents the baboon image as a superposition of $512$ simple images, each one only showing horizontal/vertical stripes.  Let us calculate a low-rank reconstruction of the image with $1$ and with $20$ components:
    baboon.1  <-  sweep(baboon.svd$u[, 1, drop=FALSE], 2, 
                  baboon.svd$d[1], "*") %*%
                       t(baboon.svd$v[, 1, drop=FALSE])
    
    baboon.20 <-  sweep(baboon.svd$u[, 1:20, drop=FALSE], 2, 
                baboon.svd$d[1:20], "*") %*%
                       t(baboon.svd$v[ , 1:20, drop=FALSE])

resulting in the following two images:

On the left we can easily see the vertical/horizontal stripes in the rank-1 image.
Let us finally look at the "residual image", the image reconstructed (as above, code not shown) from the $20$ rank-one images with the lowest singular values. Here it is:

Which is quite interesting: we see the parts of the original image that are difficult to represent as superposition of vertical/horizontal lines,  mostly diagonal nose hair and some texture, and the eyes!
