Given One Partition of a Matrix, What is the Best Way to Find a Second Partition that Ensures the Matrix is Nonsingular?

Suppose I have a matrix C, dimension m x n, m < n, and that rank(C) = m.

I wish to find a marix V, dimension n-m x n, such that the square matrix T = [C ; V] is nonsingular.

What is a good criterion to use to select a V and how can one find that V that meets the crtierion, bearing in mind that T will be used for caculations like T\A*T?

For example, suppose

>> C = [1 2 3 4 5;5 6 7 8 9];
>> rank(C)

ans =

     2

Then one possiblity is to define V by the null space of C

>> V1=null(C)'; T1 = [C;V1]; rank(T1)

ans =

     5

Is this guaranteed to work for all C from a theoretical standpoint? If so, is this appoach too limiting insofar as there are solutions for V that don't live in the null space of C? Is there a possibiltiy of problem with null(C) if C is close to not being full rank?

For this simple example, another solution can be found by inspection:

>> V2 = [zeros(3,2) eye(3)];T2 = [C;V2]; rank(T2)

ans =

     5

V2 has some appeal because it looks "simple" with only three non-zero elements that are all untiy. But can such a "simple" matrix be found when the C doesn't easily lend itself to inspection?

Neither solution looks all that appealing with respect to rcond

>> [rcond(T1) rcond(T2)]

ans =

   3.2104e-02   8.3333e-03

though maybe these aren't all that bad.

In summary, are there any suggestions on how to find V?

Answer 1

[David Goodmanson]

Hi Paul,

First of all, if C is not close to full rank, there can be numerical problems with most any calculation involving C. Ignoring that, then every solution V does have to include all of null(C), and every V that does so will have rank n. note added> And M*null(C) where M is nonsingular is also a solution as Matt has mentioned.

But V does not have to fully " live in the null space of C ", in the sense that adding linear combinations of the rows of C to any of the rows of V still gives a solution.

The " simple " V2 approach does not work because without knowing C you can never be sure that some linear combination of the rows of V2 or (some similar matrix) does not reproduce one of the rows of C. In that case the rank of the full matrix will not be n.

C = [[1 2 3 4 5]
     [5 6 7 8 9]];
V = null(C)'
A = [C; V];
cA = cond(A)
rcA = rcond(A)

  cA = 17.5328
  rcA = 0.0321

Your condition numbers are not that bad, but that has something to do with scaling. The average value of the rows of C are 3 and 8. The null space matrix V has the property that V'V = eye(n-m), so the values of V are less than 1.

C = [[1 2 3 4 5]
     [5 6 7 8 9]];
V = 5*null(C)'
then
cA = 10.8681
rcA = 0.0474

which helps a little. As opposed to

C = [100*[1 2 3 4 5]
     100*[5 6 7 8 9]];
V = null(C)'
then
  cA = 1.7533e+03
  rcA = 4.1414e-04    

so you have to make sure that the size of C and the size of V are similar.

[Paul]

Can you explain what you mean by: "every solution V does have to include all of null(C)," For example, can you show how V2 (V2' ?) includes all of null(C)?

I assume that in our second and third examples you really meant

V = null(C)'

As for the V2, I was hoping that there might be some systematic way of coming up with "simple," whatever that means, version of V.

I'm not sure what's meant by "size of C and the size of V." What does size mean in this context? Max singular value? Any suggestions on the best way to scale V1?

[David Goodmanson]

HI Paul,

[1] like most people I tend to think of matrix linear independence much more in terms of columns than rows so I hope you don’t mind if that is done in this part. So now C is 5x2, W = null(C’) is 5x3 and A = [C W]. Note that Matlab null outputs orthonormal columns, so that the null space is described by 3 perpendicular column vectors, perpendicular to the columns of C as well. W of course has rank 3.

In order for A to be rank 5, the five columns of A must be linearly independent. Suppose that V = W*M where M is 3x3 and nonsingular. Multiplication by M on the right means that each column of V is a linear combination (LC) of columns of W. The columns of V are still linearly independent but not necessarily orthogonal to each other. Still V has rank 3. So V has to equal W*M for some nonsingular M, which what I meant by including all of W.

But in addition each column of V can contain a LC of columns of C, since that does not change the rank of A. So the full V is

V = W*M  + [ for each column of V, an arbitrary LC of columns of C ]

Of course the arbitrary LC can be zero if we wish, leaving just W*M.

For V2, that LC is definitely not zero. But now for the test of V2 to see if I will produce an A of rank 5:

Since W has orthonormal columns, W’*W = I. Then

    W’*V2 = M + W’*[ … arbitrary LC of columns of C ]

But the second term is zero by the definition of the null space W.

So W’*V2 = M and M has to be nonsingular. In your case

C = [[1 2 3 4 5]
     [5 6 7 8 9]]'
W = null(C')
A = [C W]
V2 = [zeros(3,2) eye(3)]'
M = W'*V2
rank(M)
  ans = 3   

nonsingular, so it works.

*****Back to rows for C and V*****

[2] yes I did mean null(C)' and that got changed.

[3] I don't believe a simple version of V is possible without knowledge of C, for the reason I stated: you have to be sure that no linear combination of rows of V reproduces a row of C. Of course with the V2 you have, the odds are low that that would happen, and there is basically no chance at all if you chose V2 as matrix of random numbers.

[4] The 'size of V and size of C' was not a good choice of words. I meant the [ average size of the elements of rows of C] vs. [average size of the elements of rows of V ], or let's say norm(each row of V) compared to norm(each row of C). Matt has shown a good way to scale V.

[Paul]

[1] Returning to rows, it sounds like you (and Matt J.) are both saying that every V that makes T full rank will be of the form:

V = P*C + M*null(C)'   % eqn (1)

where M is non-singular. So the problem is then finding a suitbable P and M that satisfy some other criterion, like minimizing cond(T) or some (ill-defined as yet) criterion on the "simplicity" of V (all assuming that null(C) doesn't break down for C that's near rank-deficient). I think that general form of the solution is what's needed if if the criteron can be stated as an optimization over P and M.

[4] Matt J. showed how to scale T, but that's not allowed because that will also change the C partion in T. As stated in the original question the first m rows of T must be C.

[5] I will never again use C as a variable in a post. The autocorrect of (C) to a copyright symbol is driving me crazy.

[David Goodmanson]

what I called A I should have called T, sorry.

Staying with rows:

I agree with all of this.

eqn(1) is highly useful for relating conditions on T to possibly easier-to-define conditions on P and M, whichever direction that goes.

As far as cond and rond, if it matters (which it certainly doesn't in the example) you could normalize the rows of V to be the mean (arithmetic? geometric? harmonic?) of the existing norms of the rows of C. Doing one's best.

Yes, (C) sucks, Or maybe I should say that © sucks. I spend a lot of time typing (CC) and then deleting one of the Cs.

[Paul]

I've already played around with fmincon to optimize P and M for different objective functions on T, e.g. cond(T), and I at least was able to get some answers for this simple example. It's a constrained optimization because M must be nonsingular. Any thoughts on how to express that crteria as a constraint for fmincon? Fmincon only allows for nonlinear constraints of the form c(M) <= 0 or c(M) == 0. Unfortunately, c(M) < 0 isn't an option, or I would have done something like c(M) = -min(svd(M)) or something like that. I ended up using the equality constraint c(M) = rank(M) - (n - m) == 0, which worked but doesn't feel like the best approach.

[Bruno Luong]

But I did give the form of the solution

xd = smin + rand(p,1)*(smax-smin); % you could flip randomly the sign of xd if you want
M = diag(xd);

to obtain best possible COND(T). Why using FMINCON for such problem where a simple form has explicit formula?

[Paul]

I was using FMINCON because I was looking at different objective functions on T or V. Also, I misunderstood (or didn't fully understand) your answer. I did compare your solution given immediately above with the result of FMINCON. You beat FMINCON by a whisker. Here's the code, in case anyone is interested:

C = [1:5;5:9];
[m,n] = size(C);
p = n- m;
s = abs(svd(C));
smin = min(s);
smax = max(s);

V1 = null(C).'; T1 = [C;V1];
V2 = [zeros(p,m) eye(p)]; T2 = [C;V2];

rng(10);
xd = smin + rand(p,1)*(smax-smin); % you could flip randomly the sign of xd if you want
M3 = diag(xd); P3 = zeros(p,m);
V3 = P3*C + M3*V1; T3 = [C;V3];

nonlcon = @(x) (deal([],rank(x(:,n-m:n))-(n-m)));
x0 = [zeros(p,m) eye(p)];

obj4 = @(x)(cond([C;x(:,1:n-m-1)*C + x(:,n-m:n)*V1]));

x = fmincon(obj4,x0,[],[],[],[],[],[],nonlcon);
P4 = x(:,1:m); M4 = x(:,m+1:end); V4 = P4*C + M4*V1; T4 = [C;V4];

answers = {C,T1,T2,T3,T4};
Rank = cellfun(@rank,answers)
Cond = cellfun(@cond,answers)

And the results:

Rank =

     2     5     5     5     5


Cond =

     1.086814306975882e+01     1.753275524671700e+01     1.211034282587723e+02     1.086814306975881e+01     1.086814306975892e+01

Answer 2

[Matt J]

V must be of the form A*C+B*null(C).' where B is a non-singular matrix.

I don't think the rconds in your example are bad at all, but you can improve them somewhat using row normalization:

>> rcond(normalize(T1,2,'norm'))

ans =

    0.0516

>> rcond(normalize(T2,2,'norm'))

ans =

    0.0219

[Paul]

Did you mean A*null(C)' ? Note that:

>> null(C')

ans =

  2×0 empty double matrix

Also, I agree that any V of the form A*null(C)' with A rull rank will work, but it's not clear to me that those are the only solutions. After all, V2 makes T full rank, but is there an A s.t. A*V1 = V2?.

>> A=V2/V1;
>> A*V1

ans =

  -2.0000e-01  -2.0000e-01   8.0000e-01  -2.0000e-01  -2.0000e-01
  -2.7756e-17  -1.0000e-01  -2.0000e-01   7.0000e-01  -4.0000e-01
   2.0000e-01   2.7756e-17  -2.0000e-01  -4.0000e-01   4.0000e-01

I don't see how A*V1 = V2 can be the case at all insofar as the first column of V2 is zeros(3,1). Am I missing something?

I can't normalize T1 that way, because that also changes the rows of C, which is not allowed.

[Matt J]

Sorry, I meant that V must be of the form,

V=A*C+B*null(C).'

where B is a non-singular matrix. Because the row spaces of C and null(C).' are orthogonal complements, this is the only way that the rows of T will span all of , which is a requirement for non-singularity.

I can't normalize T1 that way, because that also changes the rows of C, which is not allowed.

Why not? Who is making the rules? If you can't, then you are stuck with whatever conditioning the rows of C limit you to.

Answer 3

[Bruno Luong]

C = [1 2 3 4 5;5 6 7 8 9];

[m,n] = size(C);
[Q,R] = qr(C.');
p=n-m;
X = rand(n,p);
% X(m+1:n+1:end) = X(m+1:n+1:end) + 0.5; % uncomment this to reenforce numerical rank
V = (Q*X).';
T = [C;V]

rank(T)

The important "criteria" here is X(m+k,k) ~= 0 for all k=1,...p; which is achieved by RAND() function, (rand never returns 0).

and this warranty det(T) > 0, therefore T is full rank.

To avoid "numerical" rank getting less than n, uncomment the commented line in the code, which make sure

X(m+k,k) >= 0.5 

for all k=1,...p, (or whatever striclty positive value instead of 0.5).

[Paul]

Any suggestions on how to scale V so that T rcond(T) is as close to 1 as possible?

Side questions: I was unaware that rand never returns 0. What is the source of that statement? Does rand ever return 1?

[Bruno Luong]

The condition of T is limited by the extrem singular values of C. C is given you cannot do anything to make it closer to 1, barely you can generate T that is not degrade it, like this:

[m,n] = size(C);
[Q,R] = qr(C.');
p=n-m;

s = abs(svd(C));
smin = min(s);
smax = max(s); % smax/smin is the limits of cond(T)
xd = smin + rand(p,1)*(smax-smin); % you could flip randomly the sign of xd if you want
X = [zeros(m,p); diag(xd)];

V = (Q*X).';
T = [C;V]

smax/smin
cond(T)

[Bruno Luong]

"Does rand ever return 1?"

No. RAND returns number in open interval (0,1)