# Orthogonal Projections, exp, leastsquares

## Orthogonal Projections:

Compute the orthogonal projection of the vector (2,1,-1,3) onto the plane generated by (-1,2,1,1) and (2,-1,-1,1).

Orthonormalize the basis for the plane and then project.

> u1 := vector([-1,2,1,1]); u2 := vector([2,-1,-1,1]); u3 := vector([2,1,-1,3);

```                             u1 := [ -1, 2, 1, 1 ]

u2 := [ 2, -1, -1, 1 ]

u3 := [ 2, 1, -1, 3 ]```
> v1 := evalm(u1/sqrt(7));
```                              1/2       1/2       1/2       1/2
v1 := [ - 1/7 7   , 2/7 7   , 1/7 7   , 1/7 7    ]```
> v2 := evalm( u2 - innerprod(u2,v1)*v1 );
`                        v2 := [ 10/7, 1/7, -3/7, 11/7 ]`
 Now normalize to make it of length 1.

> v2 := evalm(v2/sqrt(sum(v2[i]^2,i=1..4)));

```                  10    1/2           1/2            1/2          1/2
v2 := [ --- 231   , 1/231 231   , - 1/77 231   , 1/21 231    ]
231```
 finally project u3 onto v1 and v2 to get the answer:

> sol := innerprod(u3,v1)*v1 + innerprod(u3,v2)*v2;

```                                  1/2       19        1/2
sol := 2/7 7    v1 + ---- v2 231
77```
> sol := evalm(sol);
```                                  24
sol := [ ----, 9/11, -5/11, 3 ]
11```
 Now watch this:

> A := augment(u1,u2);

```                                     [ -1   2 ]
[        ]
[  2  -1 ]
A := [        ]
[  1  -1 ]
[        ]
[  1   1 ]```
> SOL := evalm( A &* leastsqrs(A,u3) );
```                                  24
SOL := [ ----, 9/11, -5/11, 3 ]
11```

Hey did you see that? It can be done in only one line!

## The exponential of a matrix.

### Problem:

Compute the exponential of the matrix B where:

> B := matrix(3,3,[1,0,0,0,1,2,0,2,1]);

```                                     [ 1  0  0 ]
[         ]
B := [ 0  1  2 ]
[         ]
[ 0  2  1 ]```

#### Solution:

You first need to diagonalize B. We get a diagonal matrix G and e non singular matrix P such that, B = P G P^(-1)

> eigensys := eigenvects(B);

```   eigensys :=

[-1, 1, {[ 0, -1, 1 ]}], [3, 1, {[ 0, 1, 1 ]}], [1, 1, {[ 1, 0, 0 ]}]```
> G := diag(-1,3,1);
```                                    [ -1  0  0 ]
[          ]
G := [  0  3  0 ]
[          ]
[  0  0  1 ]```
> P := augment([0,-1,1],[0,1,1],[1,0,0]);
```                                    [  0  0  1 ]
[          ]
P := [ -1  1  0 ]
[          ]
[  1  1  0 ]```
> evalm(P &* G &* inverse(P));
```                                  [ 1  0  0 ]
[         ]
[ 0  1  2 ]
[         ]
[ 0  2  1 ]```
 We got 'em.

> ExpB := evalm(P &* diag(exp(-1),exp(3),exp(1)) &* inverse(P));

```           [ exp(1)               0                           0             ]
[                                                                ]
ExpB := [    0     1/2 exp(-1) + 1/2 exp(3)   - 1/2 exp(-1) + 1/2 exp(3) ]
[                                                                ]
[    0    - 1/2 exp(-1) + 1/2 exp(3)   1/2 exp(-1) + 1/2 exp(3)  ]```
> evalf(");
```                   [ 2.718281828       0            0      ]
[                                       ]
[      0       10.22670818  9.858828739 ]
[                                       ]
[      0       9.858828739  10.22670818 ]```

Link to the commands in this file
Carlos Rodriguez <carlos@math.albany.edu>