Projection using Modified Gram-Schmidt orthogonality

Hello,

I need the Modified Gram-Schmidt orthogonalization method in my Research.

I wrote the following code for the projection using the Classic Gram-Schmidt:

function[Xp] = Project(A,B)
Xp = [] ;
u1 = B;
for i = 1:1:6
u2 = A(i,:)- (A(i,:)*u1)/(u1'*u1) * u1';
Xp = [Xp;u2] ;
end


end

I faced problems to convert the Modified Gram-Schmidt orthogonalization method into MATLAB code, which is illustrated in the following link https://en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process

under section Numerical stability.

Can anyone help me in this problem please?

[Torsten]

The question is what you want.

The two codes do very different things:

The first code orthonormalizes A, the second projects A on the subspace spanned by B.

[M]

@Torsten the second code uses the Classic Gram-Schmidt orthogonality to project A on B , my question is that I want to project A on B but using the Modified Gram-Schmidt.

I found the first code but i dont know how to apply it on the projection.

Thanks

[Torsten]

As said, projection is not normalization. To project a set of vectors A on a given vector B, you don't need Gram-Schmidt - neither classical nor modified.

You only need the projection formula

A_projected(i,:) = A(i,:)*B/(B.'*B) * B.';

[M]

@Torsten but I want it to be orthogonalized at the same time!!, Thats why I am using the Gram-Schmidt

[M]

@Torsten see in the following link :

https://en.wikipedia.org/wiki/Orthogonalization

"Methods for performing orthogonalization include:

• Gram–Schmidt process, which uses projection"

[M]

I wrote this code based on the classical Gram-Schmidt

function[Xp] = Project(A,B)
Xp = [] ;
u1 = B;
for i = 1:1:6
u2 = A(i,:)- (A(i,:)*u1)/(u1'*u1) * u1';
Xp = [Xp;u2] ;
end
end

[M]

@Torsten but there are problems in statbility in the classical Gram-Schmidt , that's why I am looking to use the modified one.

[Torsten]

As said, projecting a set of vectors (rows of A) on one other vector (B) doesn't produce any problems.

Maybe the real problem you are trying to solve is different from what you do, but my abilities as clairvoyant are limited. So you should describe what you are trying to do. If it is really projecting a set of vectors on one given vector, you don't need Gram-Schmid - you only need the projection formula supplied. If you want to project a vector on the span of a set of vectors, the problem will be different.

[M]

@Torsten Ok, in my research there is a matrix that should be orthogonalized with a certain vector ( I want to use this info in the cost function). I searched on google about the vector orthogonalization technique and I found in the following link https://en.wikipedia.org/wiki/Orthogonalization

"Methods for performing orthogonalization include:

Gram–Schmidt process, which uses projection"

And now i am trying to apply this method, it uses projection to performe orthogonalization, I tried the classic method and it gives a resonable answers but as i said i want to try the modified one , but i need a help in it!!

[Torsten]

The statement "there is a matrix that should be orthogonalized with a certain vector" makes no sense to me. Could you elaborate ?

If it means that you want to extract the portion of each row of A that is orthogonal to this certain vector, then your function "Project" does the job. But since one vector is compared with only one other vector, there is no such thing as "modified Gram Schmid" for this application.

[M]

@Torsten In my problem there is a matrix A is always orthogonalized with vector B (this is the nature of the problem), I want to take advantage from this info and apply orthognalization methods as a cost function

[M]

@Torsten but if you check the section

Numerical stability

in the following link : https://en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process

you will find something different from what you said!

[Torsten]

I already wrote that I don't understand what you mean by

In my problem there is a matrix A is always orthogonalized with vector B.

As a consequence, I also don't understand what you mean by

I want to take advantage from this info and apply orthognalization methods as a cost function

Maybe someone else in the forum gets it. Or you make an attempt to explain better.

[M]

@Torsten I just want to convert the section Numerical stability in the following link which illustrates the modified Gram-Schmidt : https://en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process

to matlab code in any example , because it faild with me.

This is all what I want, And then I will apply it to my variables.

[Torsten]

You don't have a Gram-Schmid orthogonalization if you project several vectors separately on only one other vector. And this is what your "Project" function does. So there is no such thing as "Numerical Stability" you have to care about or that would change the projection formula.

If you want to project B on the row space of A, this would be a different thing.

But I repeat myself ...

[Matt J]

@M Is there a reason you couldn't just use orth to orthogonalize A? Why must the orthogonalization be done with Modified Gram-Schmidt specifically?

[M]

@Torsten if I want to project B on the row space of A, what will be the different , in this case is there a need to use the Gram-Schmid method?. sorry but i am pretty new in this topic. Thanks

[Torsten]

in this case is there a need to use the Gram-Schmid method?

Yes. The easiest way is to orthonormalize the row space first (maybe using modified Gram-Schmid) and then use the formula under

https://math.stackexchange.com/questions/112728/how-do-i-exactly-project-a-vector-onto-a-subspace

[Torsten]

@M

The easiest way is to orthonormalize the row space first (maybe using modified Gram-Schmid) and then use the formula under

https://math.stackexchange.com/questions/112728/how-do-i-exactly-project-a-vector-onto-a-subspace

Anything you didn't understand in my answer ? Or wasn't it what you were asking for ?

[M]

Hi @Torsten, Its not what I am looking for

[Torsten]

Then I can assure you that from what you wrote, nobody will be able to understand what you are looking for.

Try to understand the problem first before looking for a solution.

[Jan]

@M: Please stop addressing specific users by messages like "Hi @xyz do you have any idea about my question please?"

Imagine what would happen, if all users do this: The most active users would receive a huge number of notifications and find less time to post answers.

Answer 1

[Matt J]

Aorth=orth(A);  %A orthogonalized

ProjB=Aorth*(Aorth.'*B); %projection of B

[M]

@Matt J I didn't ask about general orthogonality, I need the Modified Gram-Schmidt orthogonalization in my research, and still I didn't solve the problem. Thanks

[Matt J]

But you have been asked why you need Modified Gram-Schmidt specifically and haven't told us. What does it give you that orth() does not?

[M]

@Matt J Because I am working in an analytical problem and I need to discuss the difference of using several orthogonal techniques regarding my problem in the research but I have a problem in coding the Modified Gram-Schmidt and this is what i asked about in the question . Thanks

[M]

@Matt J I don't need such as a command. I need an algorithm.

[Torsten]

but I have a problem in coding the Modified Gram-Schmidt and this is what i asked about in the question .

The modified Gram-Schmidt method was perfectly coded in the program you deleted from your question.

[M]

@Torsten, That program gives the Q and R decomposition , This is not what i need

[M]

@Torsten I need the the projection algorithim of the Modified Gram-Schmidt , same as what I implemmented in the Classic Gram-Schmidt, Same as what illustrated in the section

Numerical stability

in the following link : https://en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process

[M]

@Torsten I implemented a code for it but it didn't give good answers

[Torsten]

That program gives the Q and R decomposition , This is not what i need

Q contains the orthonormal basis from the modified Gram Schmidt procedure.

http://web.mit.edu/18.06/www/Essays/gramschmidtmat.pdf

So don't delete codes that solve your problem.

[M]

@Torsten how to use the Q to obtian the Projection of A on B ? and what the input X means?

I deleted it because i thought it is not the same. Thanks

[M]

This is the code that i have delete it

function [Q,R] =  mgs(X)
    % Modified Gram-Schmidt.  [Q,R] = mgs(X);
    % G. W. Stewart, "Matrix Algorithms, Volume 1", SIAM, 1998.
    [n,p] = size(X);
    Q = zeros(n,p);
    R = zeros(p,p);
    for k = 1:p
        Q(:,k) = X(:,k);
        for i = 1:k-1
            R(i,k) = Q(:,i)'*Q(:,k);
            Q(:,k) = Q(:,k) - R(i,k)*Q(:,i);
        end
        R(k,k) = norm(Q(:,k))';
        Q(:,k) = Q(:,k)/R(k,k);
    end
end

[Matt J]

how to use the Q to obtian the Projection of A on B ?

[Q,R]=mgs(B);

ProjA=Q*Q.'*A;

[Torsten]

how to use the Q to obtian the Projection of A on B ?

ok, we go round in circles. The projection of vector A_j on B is

A_proj(j,:) = (A(j,:)*B)/(B.'*B) * B;

and a modified Gram-Schmidt method is not necessary to compute this projection.

[Matt J]

Q contains the orthonormal basis from the modified Gram Schmidt procedure...So don't delete codes that solve your problem.

Except that it would probably be more efficient to just use Matlab's qr command if they give the same result,

A=rand(5,3);

Qmgs=mgs(A)

Qmgs = 5×3

0.0841 0.5572 -0.1058 0.5426 0.1513 -0.6618 0.7485 -0.2152 0.5957 0.0052 0.7840 0.3870 0.3718 0.0754 -0.2149

[Q,R]=qr(A,0);
 Q=Q.*sign(diag(R)')

Q = 5×3

0.0841 0.5572 -0.1058 0.5426 0.1513 -0.6618 0.7485 -0.2152 0.5957 0.0052 0.7840 0.3870 0.3718 0.0754 -0.2149

function [Q,R] =  mgs(X)
    % Modified Gram-Schmidt.  [Q,R] = mgs(X);
    % G. W. Stewart, "Matrix Algorithms, Volume 1", SIAM, 1998.
    [n,p] = size(X);
    Q = zeros(n,p);
    R = zeros(p,p);
    for k = 1:p
        Q(:,k) = X(:,k);
        for i = 1:k-1
            R(i,k) = Q(:,i)'*Q(:,k);
            Q(:,k) = Q(:,k) - R(i,k)*Q(:,i);
        end
        R(k,k) = norm(Q(:,k))';
        Q(:,k) = Q(:,k)/R(k,k);
    end
end

[M]

@Torsten also I am repeating myself,my goal is to guarantee the orthogonality between A and B. Not THE PROJECTION itself! Gram Schmidt uses the projection as a tool to get the orthogonality .

For example House Householder transformation uses reflection to perform the orthogonalization and so on ....

[Torsten]

You partition A(j,:) in two parts: its projection on B ( (A(j,:)*B)/(B.'*B) * B ) and the part orthogonal to B: ( A(j,:) - (A(j,:)*B)/(B.'*B) * B ). Their sum gives A(j,:).

[M]

@Torsten But I want to apply the modified Gram–Schmidt process !!

"Methods for performing orthogonalization include:

• Gram–Schmidt process, which uses projection"

I need to study its efficiency to solve my problem. I dont need just final answer!

[Torsten]

You need Gram-Schmidt if you want to compute a projection on a set of vectors.

So if you want to project B on A, you need Gram-Schmidt.

If you want to project A on B, you need one single line of code.

[M]

@Torsten so good , I need this part "( A(j,:) - (A(j,:)*B)/(B.'*B) * B )." in my problem. This is what we call classic Gram–Schmidt process. but there is a simple modification in the modified Gram–Schmidt it should give the same answer of this part "( A(j,:) - (A(j,:)*B)/(B.'*B) * B )." .. But the coding of it faild with me.

[M]

@Torsten Ok , suppose i want to project B on A using the modified Gram–Schmidt , How can I code that?

[M]

@Torsten " if you want to compute a projection on a set of vectors.

So if you want to project B on A, you need Gram-Schmidt. "

this is exactly what i am looking for using the modified Gram–Schmidt .

[Matt J]

@M Then see above. The solution has already been presented.

[M]

@Matt J,

[Q,R]=mgs(B);

ProjA=Q*Q.'*A;

This didnt give the same answer of the classic Gram–Schmidt "( A(j,:) - (A(j,:)*B)/(B.'*B) * B )."

[M]

it should be at least close @Matt J

[Torsten]

( A(j,:) - (A(j,:)*B)/(B.'*B) * B )." in my problem. This is what we call classic Gram–Schmidt process

This is what we call "projecting a set of vectors A onto a single vector B". There is no iterative process as in the Gram-Schmidt procedure involved. And since there is no Gram-Schmidt, there is no Modified Gram-Schmidt, either.

[M]

@Torsten suppose that A is the single vector and B is the set of vectors , Then we need the Modified Gram-Schmidt

[Torsten]

Now you want to change the roles of A and B ? Do you want to fool us ?

[Q,R]=mgs(B);
ProjA=Q*Q.'*A;

should give

(A(j,:)*B)/(B.'*B) * B

not

A(j,:) - (A(j,:)*B)/(B.'*B) * B 

[M]

@Matt J , By subctracting the A from ProjA , We can get the orthogonal part , right?

[M]

@Torsten Nope, From the beginning I want any general code that do the orthogonalization by Modified Gram-Schmidt regardless my problem.

[Matt J]

The projection of column vector A onto the column space of matrix B is the unique vector P satisfying B.'*(P-A)=0. The demo below shows that the procedure from above with Q produces this point and also agrees with the more usual least squares solution for the projection B*(B\A).

B=rand(5,3);

A=rand(5,1);

[Q,R]=qr(B,0);
 Q=Q.*sign(diag(R)');
 
P1=Q*Q.'*A

P1 = 5×1

0.5921 0.7305 0.5290 0.2955 0.2067

P2=B*(B\A)

P2 = 5×1

0.5921 0.7305 0.5290 0.2955 0.2067

%Orthogonality test

B'*(A-P1)

ans = 3×1

1.0e-15 * -0.5829 -0.7772 -0.3331

[Torsten]

A = rand(5,3);
B = rand(5,1);
[Q,R]=mgs(B);
ProjA=Q*Q.'*A

ProjA = 5×3

0.0893 0.2064 0.1511 0.0225 0.0520 0.0381 0.4198 0.9709 0.7106 0.1182 0.2734 0.2001 0.0302 0.0698 0.0511

for j = 1:3
(A(:,j).'*B)/(B.'*B) * B
end

ans = 5×1

0.0893 0.0225 0.4198 0.1182 0.0302

ans = 5×1

0.2064 0.0520 0.9709 0.2734 0.0698

ans = 5×1

0.1511 0.0381 0.7106 0.2001 0.0511

function [Q,R] =  mgs(X)
% Modified Gram-Schmidt.  [Q,R] = mgs(X);
% G. W. Stewart, "Matrix Algorithms, Volume 1", SIAM, 1998.
[n,p] = size(X);
Q = zeros(n,p);
R = zeros(p,p);
for k = 1:p
    Q(:,k) = X(:,k);
    for i = 1:k-1
        R(i,k) = Q(:,i)'*Q(:,k);
        Q(:,k) = Q(:,k) - R(i,k)*Q(:,i);
    end
    R(k,k) = norm(Q(:,k))';
    Q(:,k) = Q(:,k)/R(k,k);
end
end

[M]

@Torsten Thanks. Now to get "A(j,:) - (A(j,:)*B)/(B.'*B) * B"

I have to do : A - ProjA , right?

[Matt J]

@M You must transpose A and B if you want to project rows instead of columns. Your code routine mgs() (and also matlab's qr() as we've seen) orthogonalizes the columns of the input X, not the rows.

[Torsten]

Yes, A - ProjA is the part of the columns of A orthogonal to B.

[M]

@Matt J Geart. Thank you very much

[M]

@Torsten very beneficial discussion. Thanks much.

[Torsten]

@M

I'm surprised you now found what you were searching for.

In the end, was it projecting a single vector onto a set of vectors or a set of vectors onto a single vector you were aiming at ?

[M]

@Torsten, my problem is an optimization problem. it contains equations and I need to know unknown values using several orthogonal techniques. Thanks