Program to compute partial derivatives

Question

I need to find partial derivatives for:

$L=\frac{1}{2}\sum_{i=1}^{n} w_{i}^{2}\sigma_{i}^{2} -\lambda \left( \sum_{i=1}^{n} w_{i} \bar{r_{i}} - \bar{r} \right) -\mu \left( \sum_{i=1}^{n} w_{i}-1 \right)$

with 5 variables, I get 5 partial derivatives $\frac{\partial L}{\partial w_{1}}$, $\frac{\partial L}{\partial w_{2}}$, $\frac{\partial L}{\partial w_{3}}$, $\frac{\partial L}{\partial \lambda}$ and $\frac{\partial L}{\partial \mu}$ with which I need to solve for $w_{1}$, $w_{2}$ and $w_{3}$. For all $i$, $\sigma_{i}$ is known. I don't want to manually differentiate; it could be more cumbersome with more restrictions. Later, I need to find the point where the minimum variance is located; i.e. to find out the optimal values for $\lambda$, $\mu$ and $w_{i}$ for all $i$ where $i\in\{1,2,3\}$. I need to vary the 5 variables to find the optimal solutions (a bit like in Excel but no such thing currently available) and then show the frontier graphically. So my question is which program would you use to solve such thing?

Answer 1

Just a remark: is there an inequality constraint $\sum_{i=1}^{n} w_{i} \leq 1$ in your system? I merely ask because $\mu$ is usually the KKT multiplier for inequality constraints. If indeed there is an inequality constraint, you will need to satisfy more conditions than just $\nabla L = 0$ to attain optimality (i.e. dual feasibility, complementarity slackness conditions etc.) In which case, you'd be better off using a proper optimization solver. I don't know if $\bar r_{i}, \bar r$ are constants and if the Hessian is positive semidefinite, but if so, this looks like a quadratic program (QP), and you can get solvers for this type of problem (e.g. quadprog in MATLAB).

But, if you know what you're doing.... here are some ways to get derivatives.

• Finite differencing. You did not mention if you wanted numerical or exact derivatives. If accuracy isn't too big a deal, a finite difference perturbation is the easiest option. Choose a small enough $h$ for your application and you're all set. http://en.wikipedia.org/wiki/Numerical_differentiation#Finite_difference_formulae

• Complex methods. If you want a more accurate derivative, you can also calculate a complex derivative. http://en.wikipedia.org/wiki/Numerical_differentiation#Complex_variable_methods. Essentially, the idea is this: $$F'(x) \approx \frac{\mathrm{Im}(F(x + ih))}{h} $$ Any programming language that implements a complex number type (e.g. Fortran) can be used to implement this. The derivative you get will be a real number. This is far more accurate than finite differencing, and for a sufficiently small $h$, you'll get a derivative that's pretty close to the exact derivative (to the limit of your machine precision). Note that you can only get 1st order derivatives using this method; it cannot be chained. Some people do interpolation to calculate 2nd order derivatives, but all the nice properties of the complex method are lost in that approach.

• Automatic differentiation (AD). This is the fastest and most accurate technique for obtaining numerical derivatives of any order (they can be made accurate to machine precision), but it's also the most complicated from a software point of view. Most optimization modeling languages (like AMPL or GAMS) provide AD facilities to solvers. This is a whole topic unto itself, but in short, if you're using MATLAB, you can use the INTLAB toolbox to quickly and easily calculate the derivatives you need. See this page for options: http://www.autodiff.org/?module=Tools

But for a system as small and as simple as yours, I'd just do it by hand. Or use a symbolic tool like Maxima (free) or Sage (also free, front end to Maxima).

Answer 2

Here is the example of implementation with R. R is flexible in regards of combining symbolic differentiation and numerical optimisation. This means that you can get from symbolic expression to function used in numerical optimisation quite fast. The cost is of course speed, since hand-written functions will certainly work faster.

So first create the expression:

n <- 3

wi <- paste("w",1:n,sep="")

sigma<-paste("sigma",1:n,sep="")

lfun <- paste(paste(wi,sigma,sep="^2*"),collapse="+")

rr <- paste("r",1:n,sep="")

res1 <- paste(paste(wi,rr,sep="*"),collapse="+")

res2 <- paste(wi,collapse="+")

optfun <- paste("1/2*(",lfun,")-lambda*(",res1,"-barr)-mu*(",res2,"-1)",collapse="")

vars <- c(wi,"lambda","mu")

optexpr <- parse(text=optfun)

Note that I only did string manipulation which is pretty readable. The result is:

> optexpr
expression(1/2*( w1^2*sigma1+w2^2*sigma2+w3^2*sigma3 )-lambda*( w1*r1+w2*r2+w3*r3 -barr)-mu*( w1+w2+w3 -1))
attr(,"srcfile")
<text> 

Now use symbolic differentiation to get the equations which we need to solve:

> gradexpr <- lapply(vars,function(l)D(optexpr,name=l))
> gradexpr
[[1]]
1/2 * (2 * w1 * sigma1) - lambda * r1 - mu

[[2]]
1/2 * (2 * w2 * sigma2) - lambda * r2 - mu

[[3]]
1/2 * (2 * w3 * sigma3) - lambda * r3 - mu

[[4]]
-(w1 * r1 + w2 * r2 + w3 * r3 - barr)

[[5]]
-(w1 + w2 + w3 - 1)

Now each of the variables is separate variable. For numerical optimisation it makes sense to combine indexed w, r, and sigma into vectors. This involves a little black magic, but it is not hard. We need the following function:

subvars <- function(expr,tb) {
    if(length(expr)==1) {
        nm <- deparse(expr)
        if(nm %in% tb[,1]) {
            e <- tb[tb[,1]==nm,2]
            return(parse(text=e)[[1]])
        }
        else
          return(expr)
    }
    else {
        for (i in 2:length(expr)) {
            expr[[i]] <- subvars(expr[[i]],tb)
        }          
    }       
    return(expr)
}

This function substitutes the expressions according to substitution table. Here is the end result.

> subtable <- rbind(cbind(wi,paste("w[",1:n,"]",sep="")),
                  cbind(sigma,paste("sigma[",1:n,"]",sep="")),
                  cbind(rr,paste("r[",1:n,"]",sep=""))  
                  )
> eqs <- lapply(gradexpr,function(l)subvars(l,subtable))

> eqs
[[1]]
1/2 * (2 * w[1] * sigma[1]) - lambda * r[1] - mu

[[2]]
1/2 * (2 * w[2] * sigma[2]) - lambda * r[2] - mu

[[3]]
1/2 * (2 * w[3] * sigma[3]) - lambda * r[3] - mu

[[4]]
-(w[1] * r[1] + w[2] * r[2] + w[3] * r[3] - barr)

[[5]]
-(w[1] + w[2] + w[3] - 1)

Now write a function which evaluates the equations.

eqs.fun <- function(w,sigma,r,barr,lambda,mu) {
    cl <- match.call()
    args <- as.list(cl)
    args <- args[-1]
    env <- parent.frame()
    args <- lapply(args,eval,envir=env)  
    sapply(eqs,function(l)eval(l,env=args))
}

This is a human readable format, where we supply all the neccessary data to the equations. For optimisation we need a variant of this function where all the variables and the data is separated:

nl.eqs.fun <- function(p,sigma,r,barr) {
    eqs.fun(w=p[1:3],sigma=sigma,r=r,barr=barr,lambda=p[4],mu=p[5])
}

And now we can solve the equations. For this use function nleqslv from the package nleqslv:

nleqslv(c(0,0,0,0,0),nl.eqs.fun,nl.jac.fun,sigma=c(1,1,1),r=c(1,0.5,0.5),barr=1)
$x
[1]  1.000000e+00 -1.254481e-16  2.403703e-16  2.000000e+00 -1.000000e+00

$fvec
[1] -3.330669e-16 -5.551115e-16 -1.110223e-16  0.000000e+00 -2.220446e-16

$termcd
[1] 1

$message
[1] "Function criterion near zero"

$scalex
[1] 1 1 1 1 1

$nfcnt
[1] 1

$njcnt
[1] 1

We can confirm that the solution is correct with solve.QP from the quadprog package, since our problem is quadratic optimisation problem.

solve.QP(Dmat=diag(c(1,1,1)),dvec=rep(0,3),Amat=cbind(c(1,0.5,0.5),c(1,1,1)),bvec=c(1,1),meq=2)
$solution
[1]  1.000000e+00  0.000000e+00 -1.110223e-16

$value
[1] 0.5

$unconstrained.solution
[1] 0 0 0

$iterations
[1] 3 0

$Lagrangian
[1] 2 1

$iact
[1] 1 2

In this case approach described is an overkill, since there is readily available package to solve the problem. But it is not that hard to extend it to any optimisation function. Also you can calculate the jacobian of the equations very easy.

The main advantage is that this code is much easier to debug, and if you want to change your optimisation function, you only need to change the expressions, conversion to numerical functions is automatic (if the data remains the same).

Answer 3

A bit late, if you want to use Python you have several choices:

1. You can employ Automatic Differentiation, which essentially uses the chain-rule and a look-up table to differentiate for you. Three packages I know of which do this are:

   a. OpenOpt's FuncDesigner (http://openopt.org/FuncDesigner)

   b. Theano which additionally optimizes your code (including compiling to the GPU). However, a major caveat is that in order to do its magic, it hides a lot from you (personally not my cup of tea).

   c. ScientificPython (one of its many modules)

2. You can do symbolic differentiation with SymPy (which does a large number of other things as well).

With at least 1.a, 1.c, and 2 you can get your answers and then use the answers in whatever your favorite language happens to be.