# Lecture 7, direct methods for constrained optimization

Direct methods for constrained optimization are also known as *methods of feasible directions*

## Feasible descent directions

Let $S\subset \mathbb R^n$ ($S\neq \emptyset$ closed) and $x^*\in S$. 
**Definition:** The set
$$ D = \{d\in \mathbb R^n: d\neq0,x^*+\alpha d\in S \text{ for all } \alpha\in (0,\delta) \text{ for some } \delta>0\}$$
is called the cone of feasible directions of $S$ in $x^*$.

**Definition:** The set 
$$ F = \{d\in \mathbb R^n: f(x^*+\alpha d)<f(x^*)\text{ for all } \alpha\in (0,\delta) \text{ for some } \delta>0\}$$
is called the cone of descent directions.

**Definition:** The set $F\cap D$ is called the cone of feasible descent directions.



**(Obvious) Theorem:** Consider an optimization problem 
$$
\begin{align}
\min &\  f(x)\\
\text{s.t. }&\ x\in S
\end{align}
$$
and let $x^*\in S$. Now if $x^*$ is a local minimizer **then** the set of feasible descent directions $F\cap D$ is empty.

## Idea for the methods of feasible descent directions

1. Assume a feasible solution $x$.
2. Find a feasible descent direction $d\in D\cap F$.
3. Determine the step length to the direction $d$
4. Update $x$ accordingly.




# Rosen's projected gradient method

Assume a problem with linear equality constraints
$$
\min f(x)\\
\text{s.t. }Ax=b,
$$
where $A$ is a $lxn$ matrix ($l\leq n$) and $b$ is a vector.

Let $x$ be a feasible solution to the above problem.

It holds that:

$d$ is a feasible direction *if and only if* $Ad=0$

Thus, the gradient $\nabla f(x)$ is a feasible descent direction, if 
$$ A\nabla f(x)=0.$$

This may or may not be true.

However, we can project the gradient to the set of feasible descent directions
$$ \{d\in \mathbb R^n: Ad=0\},$$
which now is a linear subspace.



### Projection

Let $a\in \mathbb R^n$ be a vector and let $L$ be a linear subspace of $\mathbb R^n$. Now, the following are equivalent
* $a^P$ is the projection of $a$ on $L$,
* $\{a^P\} = \operatorname{argmin}_{l\in L}\|a-l\|$, and
* $a^P\in A$ and $(a-a^P)^Tl=0$ for all $l\in L$.

## Projected gradient

The projection of the gradient $\nabla f(x)$ on the set $\{d\in \mathbb R^n: Ad=0\}$ is denoted by $\nabla f(x)^P$ and called the *projected gradient*. 

Now, given some conditions, the projected gradient gives us a feasible descent direction.



## How to compute the projected gradient?

There are different ways, but at this course we can use optimization. Basically, the optimization problem that we have to solve is
$$
\min \|\nabla f(x)-d\|\\
\text{s.t. }Ad=0.
$$

Since it is equivalent to minimize the square of the objective function $\sum_{i=n}\nabla_i f(x)^2+d_i^2-2\nabla_i f(x)d_i$, we can see that the problem is a quadratic problem with inequality constraints,
$$
\min \frac12 d^TId-\nabla f(x)^Td\\
\text{s.t. }Ad=0
$$
which means that we just need to solve the system of equations (see e.g., https://en.wikipedia.org/wiki/Quadratic_programming#Equality_constraints)
$$
\left[
\begin{array}{cc}
I&A^T\\
A&0
\end{array}
\right] 
\left[\begin{align}d\\\lambda\end{align}\right]
= \left[ 
\begin{array}{c}
\nabla f(x)\\
0
\end{array}
\right],
$$
where I is the identity matrix, and $\lambda$ are the KKT multipliers.

### Code in Python

#### A function for projecting a vector to a linear space defined by $Ax=0$.


```python
import numpy as np
def project_vector(A,vector):
    #convert A into a matrix
    A_matrix = np.matrix(A)
    #construct the "first row" of the matrix [[I,A^T],[A,0]]
    left_matrix_first_row = np.concatenate((np.identity(len(vector)),A_matrix.transpose()), axis=1)
    #construct the "second row" of the matrix
    left_matrix_second_row = np.concatenate((A_matrix,np.matrix(np.zeros([len(A),len(vector)+len(A)-len(A[0])]))), axis=1)
    #combine the whole matrix by combining the rows
    left_matrix = np.concatenate((left_matrix_first_row,left_matrix_second_row),axis = 0)
    #Solve the system of linear equalities from the previous page
    return np.linalg.solve(left_matrix, \
                           np.concatenate((np.matrix(vector).transpose(),\
                                           np.zeros([len(A),1])),axis=0))[:len(vector)]
```


```python
A = [[1,0,0],[0,1,0]]
gradient = [1,1,1]
project_vector(A,gradient)
```




    matrix([[ 0.],
            [ 0.],
            [ 1.]])



# Example

Let us study optimization problem
$$
\begin{align}
\min \qquad& x_1^2+x_2^2+x_3^2\\
\text{s.t.}\qquad &x_1+x_2=3\\
    &x_1+x_3=4.
\end{align}
$$
Let us project a gradient from a feasible point $x=(1,2,3)$

Now, the matrix
$$
A = \left[
\begin{array}{ccc}
1& 1 & 0\\
1& 0 & 1
\end{array}
\right]
$$.


```python
import ad
A = [[1,1,0],[1,0,1]]
gradient = ad.gh(lambda x:x[0]**2+x[1]**2+x[2]**2)[0]([1,2,3])
d = project_vector(A,[-i for i in gradient])
print d
```

    [[ 2.66666667]
     [-2.66666667]
     [-2.66666667]]


### d is a feasible direction


```python
np.matrix(A)*d
```




    matrix([[  0.00000000e+00],
            [  4.44089210e-16]])



### d is a descent direction


```python
def f(x):
    return x[0]**2+x[1]**2+x[2]**2
alpha = 0.001
print "Value of f at [1,2,3] is "+str(f([1,2,3]))
x_mod= np.array([1,2,3])+alpha*np.array(d).transpose()[0]
print "Value of f at [1,2,3] +alpha*d is "+str(f(x_mod))
print "Gradient dot product direction (i.e., directional derivative) is " \
+ str(np.matrix(ad.gh(f)[0]([1,2,3])).dot(np.array(d)))
```

    Value of f at [1,2,3] is 14
    Value of f at [1,2,3] +alpha*d is 13.978688
    Gradient dot product direction (i.e., directional derivative) is [[-21.33333333]]


## Finally, the algorithm of the projected gradient


```python
import numpy as np
import ad
def projected_gradient_method(f,A,start,step,precision):
    f_old = float('Inf')
    x = np.array(start)
    steps = []
    f_new = f(x)
    while abs(f_old-f_new)>precision:
        f_old = f_new
        gradient = ad.gh(f)[0](x)
        grad_proj = project_vector(A,[-i for i in gradient])#The only changes to steepest..
        grad_proj = np.array(grad_proj.transpose())[0] #... descent are here!
#        import pdb; pdb.set_trace()
        x = x+grad_proj*step
        f_new = f(x)
        steps.append(list(x))
    return x,f_new,steps
```


```python
f = lambda x:x[0]**2+x[1]**2+x[2]**2
A = [[1,1,0],[1,0,1]]
start = [1,2,3]
(x,f_val,steps) = projected_gradient_method(f,A,start,0.6,0.000001)
```


```python
print x
print f(x)
print f([1,2,3])
print np.matrix(A)*np.matrix(x).transpose()
```

    [ 2.333248  0.666752  1.666752]
    8.66666668851
    14
    [[ 3.]
     [ 4.]]


8008c88efa0276f5e136195ca53db243608bb6b1
ipynb
Jupyter Notebook
Lecture 7, Direct methods for constrained optimization.ipynb
maeehart/TIES483
cce5c779aeb0ade5f959a2ed5cca982be5cf2316
CC-BY-3.0
2019-04-26T12:46:14.000Z
2021-11-23T03:38:59.000Z
Lecture 7, Direct methods for constrained optimization.ipynb
maeehart/TIES483
cce5c779aeb0ade5f959a2ed5cca982be5cf2316
CC-BY-3.0
Lecture 7, Direct methods for constrained optimization.ipynb
maeehart/TIES483
cce5c779aeb0ade5f959a2ed5cca982be5cf2316
CC-BY-3.0
2016-01-08T16:28:11.000Z
2021-04-10T05:18:10.000Z
Qwen/Qwen-72B
1. YES
2. YES