How should the L matrix in LU decomposition be padded when m > n?

Question

I understand that L needs to be square, but all the examples that I could find for doing decomposition on non-square matrices is when the number of columns is greater than the number of rows. How should L be squared when the number of rows is greater than the number of columns?

The way Cuda's LU decomposition function works is that it overwrites the matrix A by its factors. The causes both L and U to have the dimensions of the original matrix and I need figure out how to square the L factor when A is non-square.

As a guess I am trying to set the diagonals of the L matrix returned to me by the Cuda library's LU decomposition function to the identity and the rest to zeroes, but it does not seem to be the correct choice.

Answer 1

Going by the online calculator it seems that my initial guess of padding the off diagonal elements of L with zeroes and the diagonal with ones is correct. The reason why this is not working for me is because for some reason the Cuda library rather than solving for L, just leaves the elements of the original matrix where they were. This was due to me passing arguments in the incorrect order.

Answer 2

There is an implementation at Matlab's file exchange of complete pivoting. The code looks like this. There is a part with a conditional statement towards the bottom handling non-square matrices.

function [L U p q] = lucp(A,tol,pm_opt)
%LUCP     LU factorization with complete pivoting.
%
% To compute the LU factorization under default settings:
%
%  [L U p q] = lucp(A)
%
% This produces a factorization such that L*U = A(p,q).  Vectors p and q
% permute the rows and columns, respectively.
%
% The pivot tolerance can be controlled:
%
%  [L U p q] = lucp(A,tol)
%
% The algorithm will terminate if the absolute value of the pivot is less
% than tol.
%
% Permutation matrices can be generated:
%
%  [L U P Q] = lucp(A,tol,'matrix')
%  [L U P Q] = lucp(A,tol,'sparse')
%
% The first will generate full permutation matrices P and Q such that 
% L*U = P*A*Q.  The second generates sparse P and Q.
%
% If A is sparse, L and U will be sparse.
%
% This function works on non-square matrices.
%
% Input:
%  A = matrix
%  tol = pivot tolerance (defualt is 1e-10)
%  pm_opt = permutation output options
%         = 'vector' for permutation vectors, L*U = A(p,q), defualt
%         = 'matrix' for full permutation matrices, L*U = P*A*Q
%         = 'sparse' for sparse permutation matrices, L*U = P*A*Q
%
% Output:
%  L = lower triangular factor
%  U = upper triangular factor
%  p = row permutation
%  q = column permutation
%
% Reference:
%  Algorithm 3.4.2, Matrix Computations, Golub and Van Loan. (3rd ed)
%  
% Other Implementations:
%  Gaussian Elimination using Complete Pivoting.  Alvaro Moraes.
%  http://www.mathworks.com/matlabcentral/fileexchange/25758
%
%  Gauss elimination with complete pivoting.  Nickolas Cheilakos.
%  http://www.mathworks.com/matlabcentral/fileexchange/13451
%  (Does not work with rectangular A)
%
%  Rank Revealing Code.  Leslie Foster.
%  http://www.math.sjsu.edu/~foster/rankrevealingcode.html
%  (Uses mex libraries for computation)
%


%
% 2010-03-28 (nwh) first version.
% 2010-04-14 (nwh) added more references.
% 2010-04-24 (nwh) perform final column swap so U is well conditioned
%

%
% License: see license.txt.
%

% handle optional inputs
if nargin < 2 || isempty(tol)
  tol = 1e-10;
end

if nargin < 3 || isempty(pm_opt)
  pm_opt = 'vector';
end

if strcmp(pm_opt,'vector')
  pm_flag = false;
  sp_flag = false;
elseif strcmp(pm_opt,'matrix')
  pm_flag = true;
  sp_flag = false;
elseif strcmp(pm_opt,'sparse')
  pm_flag = true;
  sp_flag = true;
else
  error('lucp:invalidOption','''%s'' is an un recognized option.',pm_opt)
end

[n m] = size(A);

% pivot vectors
p = (1:n)';
q = (1:m)';

% temp storage
rt = zeros(m,1); % row temp
ct = zeros(n,1); % col temp
t = 0; % scalar temp

for k = 1:min(n-1,m)
  % determine pivot
  [cv ri] = max(abs(A(k:n,k:m)));
  [rv ci] = max(cv);
  rp = ri(ci)+k-1;
  cp = ci+k-1;

  % swap row
  t = p(k);
  p(k) = p(rp);
  p(rp) = t;
  rt = A(k,:);
  A(k,:) = A(rp,:);
  A(rp,:) = rt;

  % swap col
  t = q(k);
  q(k) = q(cp);
  q(cp) = t;
  ct = A(:,k);
  A(:,k) = A(:,cp);
  A(:,cp) = ct;

  if abs(A(k,k)) >= tol
    rows = (k+1):n;
    cols = (k+1):m;
    A(rows,k) = A(rows,k)/A(k,k);
    A(rows,cols) = A(rows,cols)-A(rows,k)*A(k,cols);
  else
    % stop factorization if pivot is too small
    break
  end
end

% final column swap if m > n
if m > n
  % determine col pivot
  [cv ci] = max(abs(A(n,n:m)));
  cp = ci+n-1;

  % swap col
  t = q(n);
  q(n) = q(cp);
  q(cp) = t;
  ct = A(:,n);
  A(:,n) = A(:,cp);
  A(:,cp) = ct;
end

% produce L and U matrices
% these are sparse if L and U are sparse
l = min(n,m);
L = tril(A(1:n,1:l),-1) + speye(n,l);
U = triu(A(1:l,1:m));

% produce sparse permutation matrices if desired
if pm_flag
  p = sparse(1:n,p,1);
  q = sparse(q,1:m,1);
end

% produce full permutation matrices if desired
if ~sp_flag
  p = full(p);
  q = full(q);
end

end

Running this on a non-square matrix

A = rand(3,4)
tol =1e-8
[L U P Q] = lucp(A,tol,'matrix')


L =

    1.0000         0         0
    0.3337    1.0000         0
    0.7312    0.8581    1.0000


U =

    0.9502    0.1869    0.3816    0.6463
         0    0.7328    0.3114    0.2299
         0         0   -0.5118   -0.1801


P =

     0     0     1
     0     1     0
     1     0     0


Q =

     1     0     0     0
     0     0     1     0
     0     1     0     0
     0     0     0     1