MATLAB Answers

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

74 views (last 30 days)
Paul
Paul on 31 Dec 2020
Edited: Paul on 2 Jan 2021
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?

  0 Comments

Sign in to comment.

Accepted Answer

David Goodmanson
David Goodmanson on 1 Jan 2021
Edited: David Goodmanson on 1 Jan 2021
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.

  7 Comments

Show 4 older comments
Paul
Paul on 2 Jan 2021
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
Bruno Luong on 2 Jan 2021
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
Paul on 2 Jan 2021
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

Sign in to comment.

More Answers (2)

Matt J
Matt J on 31 Dec 2020
Edited: Matt J on 1 Jan 2021
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

  2 Comments

Paul
Paul on 1 Jan 2021
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
Matt J on 1 Jan 2021
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.

Sign in to comment.


Bruno Luong
Bruno Luong on 1 Jan 2021
Edited: Bruno Luong on 1 Jan 2021
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).

  3 Comments

Paul
Paul on 1 Jan 2021
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
Bruno Luong on 1 Jan 2021
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)

Sign in to comment.

Products


Release

R2019a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by