# Program to compute partial derivatives

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?

• Because $\lambda$ and $\mu$ are Lagrange multipliers, it's trivial to find derivatives with respect to them. How do you anticipate the functional form of $L$ changing so that the derivatives with respect to the $w_i$ -- which currently are dead simple to compute -- will become complicated? // Anyway, any symbolic mathematics program should be able to do the job well: MatLab, Mathematica, etc. – whuber Mar 13 '11 at 18:14
• @whuber: as far as I know, MATLAB is not symbolic. Perhaps you meant Maple? – posdef Mar 14 '11 at 10:00
• @posdef Matlab performs symbolic differentiation: mathworks.com/help/toolbox/symbolic/brvfu8o-1.html#brvfxct – whuber Mar 14 '11 at 14:59
• @whuber: cool, I didn't know that. thanks :) – posdef Mar 14 '11 at 15:04
• @posdef And there's sympy for Python. – chl Mar 14 '11 at 15:41

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).

• yes, the restriction $\sum_{i=1}^{n}w_{i}= 1$. – hhh Mar 14 '11 at 8:36

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))
[]
1/2 * (2 * w1 * sigma1) - lambda * r1 - mu

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

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

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

[]
-(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)[])
}
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
[]
1/2 * (2 * w * sigma) - lambda * r - mu

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

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

[]
-(w * r + w * r + w * r - barr)

[]
-(w + w + w - 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,mu=p)
}


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.000000e+00 -1.254481e-16 2.403703e-16 2.000000e+00 -1.000000e+00$fvec
 -3.330669e-16 -5.551115e-16 -1.110223e-16  0.000000e+00 -2.220446e-16

$termcd  1$message
 "Function criterion near zero"

$scalex  1 1 1 1 1$nfcnt
 1

$njcnt  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.000000e+00  0.000000e+00 -1.110223e-16

$value  0.5$unconstrained.solution
 0 0 0

$iterations  3 0$Lagrangian
 2 1

$iact  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). • @mpiktas: which readily available package you mean to solve the problem? – hhh Mar 14 '11 at 9:40 • @hhh, quadprog, since the problem stated is quadratic optimisation problem. – mpiktas Mar 14 '11 at 9:41 • @mpiktas It's nice to see how R can do the symbolic derivatives. But seeing how it takes an expert like yourself 2-3 pages of code to set up and test the operations, and seeing how unreadable it still is (relative to a truly symbolic environment like Mathematica) is a convincing demonstration of why one would not prefer R for such purposes. – whuber Mar 14 '11 at 15:01 • @whuber, I sacrificed a bit of readability to the scalability. Function$L$has 5 parameters and 7 constants. If we increase$n$to 10 we will have 12 parameters and 21 constants, but my fitting function will still require 4 parameters, since constants are separated to two vectors and 1 constant. I fail to see how using symbolic environment could be advantageous for$n=10$, since typing 33 individual variables and keeping track of them should be tedious. However my experience with symbolic environments is very limited, so I would be glad if I am proved wrong. – mpiktas Mar 14 '11 at 20:39 • @mpiktas In Mathematica, for example:$\text{With}\left[\{n=3\}, D\left[\frac{1}{2} \sum _{i=1}^n w_i^2 \sigma _i^2-\lambda \left(-r+\sum _{i=1}^n r_i w_i\right)-\mu \left(-1+\sum _{i=1}^n w_i\right), \left\{\text{Join}\left[\{\lambda , \mu \}, \text{Array}\left[w_{\#}\&,n\right]\right]\right\}\right]\right]$. The output is$\left( \begin{array}{c} r-r_1 w_1-r_2 w_2-r_3 w_3 \\ 1-w_1-w_2-w_3 \\ -\mu -\lambda r_1+w_1 \sigma _1^2 \\ -\mu -\lambda r_2+w_2 \sigma _2^2 \\ -\mu -\lambda r_3+w_3 \sigma _3^2 \end{array} \right)\$ – whuber Mar 14 '11 at 20:49

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.