%% Solving Commodious Linear Systems
% This is about linear systems with fewer equations than variables; A*x = b where the m -by- n matrix A has fewer rows that columns, so m < n . I have always called such systems wide or fat, but this is not respectful. So I consulted the Merriam-Webster Thesaurus and found commodious.
% This is also one my first posts to use the Live Editor™. Here is a a link to the .mlx file.
%% Commodious linear systems
% Commodious linear systems with fewer rows than columns do not have unique solutions.  It is hard to compute something if you don't quite know what it is you're trying to compute.  You can pin down the solution by adding an additional condition: minimize its Euclidean or \ell_2 norm.   This is the Moore-Penrose solution.  It is unique and has many other desirable properties.
% But recent years have seen many situations where a different condition is desirable: minimize the number of nonzero components of the solution.  Such a solution is called sparse.  Actually minimizing the number of nonzero components is combinatorically difficult.  (It is NP-hard.)  However, the MATLAB backslash operator produces a solution with at most m nonzero components.  This is often the minimum.
%% Generate examples
% Pick the dimensions.

   m = 3
   n = 10

%

   assert(m < n)
   format short

% Generate an m-by-n matrix with random one-digit positive and negative integer entries. Every time this script is executed, a different A is created.

   A = randi(18,m,n) - 9

% Generate a random right hand side with the same number of rows.

   b = randi(18,m,1) - 9

% Save a copy of the rhs.

   b_copy = b; 

%% Backslash
%
% Use backslash to compute one of many possible solutions to A*x = b.

   x = A\b 

% This x is sparse; it has only m nonzero components.
%
% How was x computed?   Compute the QR decomposition with column pivoting.  That is A*P = Q*R with orthogonal Q, right trapezoidal R, and permutation P.

   [Q,R,P] = qr(A)   

% Keep only the first m columns of R.

   R = R(:,1:m)

% Apply Q to right hand side.

   b = Q'*b   

% Back substitution solves the square, triangular system R*x = b.

   x = R\b  

% Tack some zeros on the end of x, matching the number of columns deleted from R.

   x(m+1:n) = 0

% Undo the permutations.

   x = P*x

%% Pseudoinverse
% Restore the rhs.

   b = b_copy;

% For comparison, the minimum 2-norm solution can be obtained with the Moore-Penrose pseudoinverse.  This is usually completely full; none of its components are zero.

   y = pinv(A)*b   

% Both x and y are solutions to the original system.

   Ax = A*x
   Ay = A*y

% However norm(y,2) is less than norm(x,2).

   normyx = [norm(y) norm(x)]

%
%% SVD
% To compute y without computing the pseudoinverse, use the SVD.

   [U,S,V] = svd(A)   

% Then, with appropriate parentheses, and with a backslash involving a diagonal matrix, finding the minimum norm solution is a one-liner, .

   y = V*(S\(U'*b))

%% Questions
% There are lots of things that I don't know about commodious systems.
%% 1.
%     Is there a  bound on norm(x)/norm(y)?
%% 2.
%     The condition number of A with respect to any norm could be defined by
%         cond(A) = max(norm(A*x)/norm(x))/min(norm(A*x)/norm(x))
% We know that for the 2-norm this is the ratio of singular values.  What about the 1-norm or the inf-norm?  How can we compute this?
%% 3.
%     What can be proved about the residual and the error of the solution computed by backslash?  I suspect the residual is small, but what is the error when the solution is not unique?
%% 4.
%     What is this distribution?

   kmax = 10000; 
   kappa = zeros(kmax,1); 
   for k = 1:kmax
       A = randi(18,m,n) - 9; 
       kappa(k) = cond(A);
   end 
   histogram(kappa)
   title('cond(A)') 
   ax = axis;
   text(0.7*ax(2),0.8*ax(4),sprintf('m =  %2d, n = %2d',m,n))