Matrix column updation via optimization

I have an initial matrix B of size (n x n). I have to update each columns of the matrix except the first column such that the updated matrix is orthogonal (BB^T =I) and it also satisfies the constraint (say Bx=c). Is there any existing optimization algorithm to solve it?

[KALYAN ACHARJYA]

Any sample example?

[John D'Errico]

I'm sorry, but this makes no sense.

You are asking to modify the matrix B, reducing it to an orthogonal matrix. But then your "constraint" introduces a completely new variable, x. What is x, and what are you then hoping to optimize? Such an orthogonality reduction for B is a simple one (gram-schmidt is an example), that needs nothing more than linear algebra.

My guess is you wanted to speicify some sort of constraint on the reduction to an orthogonal matrix, and you just wrote down the first thing that came to mind, even if it made no sense in context. Are x and c known vectors?

I can see other problems too. If you will modify the columns of b, not touching the FIRST column, then it is required that the first column of B must already be unit normalized. And that would say you may need to modify even the first column of B. Otherise it is impossible for B'*B to yield the identity.

[Veena Narayanan]

@KALYAN ACHARJYA@John D'Errico Thank you for the reply. I shall elaborate the question.

The vectors x and c are known, and c is a sparse vector. The initial matrix B is to be updated so as to satify the following conditions:

1) The first column of the matrix B remains unaltered (The first column of B is unit normalized).

2) The remainining (n-1) columns of the matrix B should be orthogonal to each other as well as to the first column of the matrix. (The first column needs to be fixed as i have to use it in another stage of the algorithm).

3) The matrix should be updated such that norm(Bx-c) is minimum.

4) The other (n-1) updated columns of B should also be unit norm vectors.

[Torsten]

You wrote down the optimization problem you have to solve:

min: (norm(B*x - c))^2

under the constraints

B*B^T = I
B(:,1) = given from original matrix B (should be a vector of norm 1)

I don't know if there is a theoretical solution.

Try "fmincon" to solve.

Do x and c have the same vector norm ?

[Veena Narayanan]

@Torsten I have tried fmincon but it did'nt work. Yes, the vectors x and c have the same vector norm.

[Torsten]

What do you mean by "didn't work" ?

Can you include a test code ?

[Veena Narayanan]

@Torsten Initially I tried this code by including only the orthogonality condition to check if the updated matrix B_hat satisfies the orthogonality condition. But the matrix B_hat was same as the initial matrix B_in and there was no updation.

The initial matrix B_in considered was a random matrix of size (64x64). c and x are vectors of size (64x1).

A = []; b = []; Aeq = []; beq = []; lb = []; ub = []; x0=B_in;

fun=@(B)norm(c-B*x)^2;

nonlcon = @orthon;

[B_hat,fval,exitflag,output] = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon)

---------------------------------------------------------

function [c,ceq] = orthon(B)

c = [];

ceq = norm(eye(64)- B*B','fro')^2;

----------------------------------------------------------

[Torsten]

If this is the problem you are trying to solve, it has a simple solution:

If you define

u = (x-c)/norm(x-c)

and

B = eye(64) - 2*u*u'

then B*x = c and B is an orthogonal matrix.

The difficulty arises because you want to keep the first column of the matrix B unchanged.

[Veena Narayanan]

@Torsten Thank you. But the orthogonality constraint should be satifiied retaining the first column same as that of the matrix B.

[Torsten]

But the orthogonality constraint should be satifiied retaining the first column same as that of the matrix B.

I don't understand. If you construct a matrix B' as above, it is orthogonal and satisfies B'*x = c.

But it is not related to your initial matrix B. Especially if you replace the first column of B' with the first column of B, B'*x = w and B' orthogonal will no longer be true in general.

I posted this "solution" because you don't use in your code that you want to prescribe the first column of B'. So the

B' = eye(64) - 2*u*u'

from above could be one solution from your code if run successfully.

[Veena Narayanan]

@Torsten My aim is to obtain an updated matrix B of size (n x n) that satisfies the condition Bx=c, which is orthogonal and have its first column same as that of the initial matrix. Therefore, the updation should be done such that the remaining (n-1) columns (except the first column) are mutually orthogonal as well as orthogonal to the first column.

The first column should be same as the initial matrix because i have to use that column at a later stage in my algorithm. I understand that the code above satisfies the orthogonality condition. But it cannot be used here since the first column also gets altered in the process.

[Torsten]

I know this, but your MATLAB code calling "fmincon" does not fix the first column. So I thought you relaxed the condition on B.

But see Bruno Luong's code below which seems to solve your problem.

[Veena Narayanan]

@Torsten Thanks a lot.

Answer 1

[Bruno Luong]

Using Housholder transformation, B is uniquely determined only for n=3.

I claim that my code solve the problem of

argmin(norm(B*x-c))

under constraints

B(;,1) = B1, and

B'*B = I

where B1, x, c are known.

PS: I assumeall variables are reals. For complex some transpose operations would be modified.

% Generate some fake data
n = 3;
[B,~] = qr(randn(n))

B = 3×3

-0.3607 -0.8710 -0.3336 0.4030 0.1770 -0.8979 0.8411 -0.4583 0.2872

x = randn(n,1);
c = B*x;
B1 = B(:,1);
clear B

% Try to recover B(:,2:end) (with sign of columns that is random) from B(:,1), x, and c
N = null(B1');
cp = c-B1*x(1);
xp = x(2:n);
cp = (cp'*N)';
vp = (cp-xp); % EDIT sign corrected
vp = vp/norm(vp);
P = eye(n-1)-2*vp*vp';
H = N*P;
B = [B1, H]

B = 3×3

-0.3607 -0.8710 -0.3336 0.4030 0.1770 -0.8979 0.8411 -0.4583 0.2872

[Veena Narayanan]

@Bruno Luong This code given an orthogonal matrix and also retains the first column. But the condition Bx=c is not satisfied with the updated matrix. The error between the original c and the c^ obtained using updated matrix B is large.

[Bruno Luong]

Your claims does mean anything to me.

The solution of my code meets B*x = c in the least square sense.

Since B is orthogonal |B*x|must be equal to |x| so if you provide x and c that does not have equal norm it can NEVER meet exactly the condition.

I bet if you could find an example where where |B*x-c| with B given by my code is larger than some other B with the same constraints. Until you don't you cannot show my code doesn't work.

[Torsten]

But your original B satisfied B*x = c exactly and was orthogonal.

So if you had recovered the original B, both conditions (B*B.' = I and B*x = c exactly) would hold.

In the case n=3, it's only the sign of P which is wrong.

I don't know if experimenting with the sign of P suffices for n > 3.

[Bruno Luong]

Test with n =5, the norm of B*x-c is pratically 0

% Generate some fake data
n = 5;
[B,~] = qr(randn(n))

B = 5×5

-0.4189 0.3461 -0.1046 0.7862 -0.2750 0.3966 -0.2379 0.7354 0.4692 0.1584 -0.0667 -0.8680 -0.4058 0.2722 -0.0585 -0.8141 -0.2212 0.4420 -0.1957 0.2339 0.0093 0.1459 -0.2970 0.2220 0.9171

x = randn(n,1);
c = B*x;
B1 = B(:,1);
clear B
% Try to recover B(:,2:end) (with sign of columns that is random) from B(:,1), x, and c
N = null(B1');
cp = c-B1*x(1);
xp = x(2:n);
cp = (cp'*N)';
vp = (cp-xp);
vp = vp/norm(vp);
P = eye(n-1)-2*vp*vp';
H = N*P;
B = [B1, H]

B = 5×5

-0.4189 -0.5396 0.3613 0.6277 -0.0940 0.3966 0.5942 0.1535 0.6818 -0.0351 -0.0667 0.2291 0.9007 -0.3624 0.0237 -0.8141 0.5477 -0.1846 0.0399 0.0407 0.0093 -0.0579 0.0257 0.0905 0.9938

norm(B*x-c)/norm(c)

ans = 1.9173e-16

[Torsten]

My respect - well done indeed.

[Veena Narayanan]

@Bruno Luong Thank you for the clarification. In the code given above, both c and x are of equal norm. Also |B*x|= |x|. Even then |B*x-c| gives a value approximately equal to 1. So do you mean to say, even if i consider some other B, i will not be able to obtain a smaller error value for |B*x-c|?

[Torsten]

Give us an example where the code does not work in your opinion.

In the examples I tested, it was successful and gave B with B*x = c exactly.

Of course, |x| = |c| must hold since an orthogonal matrix preserves lengths.

[Bruno Luong]

I fix one bug in this line

vp = (cp-xp);

so some conclusions must be revised.

[Veena Narayanan]

@Bruno Luong @Torsten Thank you so much. Now it works.

[Veena Narayanan]

@Bruno Luong Sir, in the code above, can you clarify why we consider the Householder matrix P associated with vp? Is it to obtain a set of orthogonal vectors?

[Bruno Luong]

The purpose of P (Householder reflection) is to map the projection of x (xp) to the projection of c. The projection is on the subspace orthogonal to span<B1> whith is span(N). The Housholder are reflected vectors mirror to span(vp); computed so that xp is mapped to cp. P is orthognal so N*P.