Solve indeterminate equation system in R

Question

I have a equation system and I want to solve it using numerical methods. I want to get a close solution given a starting seed. Let me explain.

I have a vector of constants ,X, of values:

X <- (c(1,-2,3,4))

and a vector W of weights:

W <- (c(0.25,0.25,0.25,0.25))

I want that the sum of the components of W will be (sum(W)=1), and the sum of the multiplication of X and W element by element will be a given number N (sum(W*X)=N).

Is there a easy way to do this in R? I have it in Excel, using Solver, but I need to automatize it.

Answer 1

Here is your constant and your target value:

x <- c(1, -2, 3, 4)
n <- 10

You need a function to minimize. The first line contains each of your conditions, and the second line provides a measure of how to combine the errors into a single score. You may want to change the second line. For example, you could make one error term be more heavily weighted than the other using sum(c(1, 5) * errs ^ 2).

fn <- function(w)
{
  errs <- c(sum(w) - 1, sum(x * w) - n) 
  sum(errs ^ 2)
}

The simplest thing is to start with all the weights the same value.

init_w <- rep.int(1 / length(x), length(x))

Use optim to optimize.

optim(init_w, fn)
## $par
## [1]  0.1204827 -1.2438883  1.1023338  1.0212406
## 
## $value
## [1] 7.807847e-08
## 
## $counts
## function gradient 
##      111       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

The par element contains your weights.

Answer 2

There is no unique solution for this problem. If you try other initial values for w you will most likely get different results from optim.

The problem can be formulated as solving an underdetermined system of linear equations.

A <- matrix(c(rep(1,4),x), nrow=2,byrow=TRUE)
b <- matrix(c(1,n), nrow=2)

We seek a solution that satisfies A %*% w = b but which one? Minimum norm solution? Or maybe some other one? There are infinitely many solutions. Solutions can be given using the pseudo-inverse of the matrix A. Use package MASS for this.

library(MASS)
Ag <- ginv(A)

The minimum norm solution is

wmnorm <- Ag %*% b

And check with A %*% wmnorm - b and fn(wmnorm). See the Wikipedia page System of linear equations the section Matrix solutions.

The solutions are given by

Az <- diag(nrow=nrow(Ag)) - Ag %*% A
w <- wmnorm + Az %*% z

where z is an arbitrary vector of ncol(Az) elements. And now generate some solutions and check

xb <- wmnorm
z <- runif(4)
wsol.2 <- xb + Az %*% z
wsol.2
A %*% wsol.2 - b
fn(wsol.2)

z <- runif(4)
wsol.3 <- xb + Az %*% z
wsol.3
A %*% wsol.2 - b
fn(wsol.3)

And you'll see that these two solutions are valid solutions when given as argument to fn. And are quite different from the solution found by optim. You could test this by choosing a different starting point init_w for example by init_w1 <- runif(4)/4.