Decomposition aproach to solve projection problem

gistfile1.matlab

% X         - source matrix 
% splitFunc - function handler for some matrix decomposition method, accepts X and initial vector,
%   which formed as convex combinations of vectors from X

function z = SimProPart(X, splitFunc)
    do % TODO: find a better a way to force non empty decomposition
        V = rand(1, columns(X));
        V = V / norm(V);
        [X1, X2] = splitFunc(X, sum(X * diag(V), 2));
    until (columns(X1) > 0 && columns(X2) > 0)

    x1 = SimPro(X1, 1e5, 1.e-8, [-1 -1], [], []); % TODO: don't repeat parameters lists
    x2 = SimPro(X2, 1e5, 1.e-8, [-1 -1], [], []);

    z = SimPro([x1 x2], 1e5, 1.e-8, [-1 -1], [], []);

    while any((dot(X, repmat(z, 1, columns(X))) - dot(z, z)) < 0)
        x1 = SimPro([X1 z], 1e5, 1.e-8, [-1 -1], [], []);
        x2 = SimPro([X2 z], 1e5, 1.e-8, [-1 -1], [], []);

        z = SimPro([x1 x2], 1e5, 1.e-8, [-1 -1], [], []);
    endwhile
endfunction

% splitFunc example
function [B, C] = splitMatrix(A, p)
    criteria = dot(A, repmat(p, 1, columns(A))) - dot(p, p);
    B = A(:,find(criteria <= 0));
    C = A(:,find(criteria > 0));
endfunction

% Bad case matrix
A = [
   3.3357e+01  -1.6679e+01  -1.6679e+01   1.5850e-03;
  -5.1799e+00   1.0360e+01  -5.1799e+00   2.3623e-03;
   1.0000e-02   1.0000e-02   1.0000e-02   2.0000e-02];