--- /dev/null
+function Ak = rank_k_approximation(A,k)
+ %
+ % Compute the rank-k approximation to A from the ``k`` largest
+ % singular values of A.
+ %
+ % INPUT:
+ %
+ % - ``A`` -- The matrix to approximate.
+ %
+ % - ``k`` -- The rank of the resulting matrix; i.e. the number
+ % of (large) singular values to keep.
+ %
+ % OUTPUT:
+ %
+ % - ``Ak`` -- The rank-k approximation to ``A``.
+ %
+ [U, S, V] = svds(A, k);
+ Ak = U*S*V';
+end
--- /dev/null
+U = [4/5, 3/5, 0 ; ...
+ 0 , 4/5, 3/5; ...
+ 3/5, 0 , 4/5 ];
+
+S = [ 1, 0, 0;
+ 0, 2, 0;
+ 0, 0, 3 ];
+
+V = [ 0, 1, 0;
+ 1, 0, 0;
+ 0, 0, 1 ];
+
+A = U*S*V';
+
+expected = A;
+actual = rank_k_approximation(A, 3);
+
+unit_test_equals("Full rank approximation of a matrix is itself", ...
+ expected, ...
+ actual);