Optimal arrangement of a small vector within a lager vector

Hi all,this is about the 3rd time I am posting this question as I do not have an answer yet. I simplified the original problem.

I have 2 vectors:

 a = [1 2 3 2 1];
 b = zeros(1, 11);

I want to find the best placements, 'p' (maximum 3 or 4x) and weightings, 'w', of vector 'a' inside vector 'b' so that the values of vector 'b' have the least standard deviation. All elements of 'b' must be non-zero.

As an example, a possible conceptual solution could be:

 b(p1)  = b(p1) + a * w1;  %w1 = 3;   p1 = 1:5;
 b(p2)  = b(p2  + a * w2;  %w2 = 2;   p2 = 4:8;
 b(p3)  = b(p3) + a * w3;  %w3 = 2;   p3 = 6:10;
 b(p4)  = b(p4) + a * w4;  %w4 = 3;   p4 = 9:11;

My problem is how to find the actual values of w1:w4 and p1:p4. Your help and suggestions will be very much appreciated.

[Wolfgang Schwanghart]

The standard deviation of the vector b (std(b)) won't be effected by the placement of a by the index vector p. The least standard deviation will be obtained by w = 0. I think you need to formulate your problem differently.

[Charles]

Hi Wolfgang, at w = 0, the elements of b would be 0. I specified that b should have non-zero elements. The placements and weighting of a inside b would clearly affect the standard deviation.

[Image Analyst]

I'm with Wolfgang. I can't see how the placement of a within b affects the standard deviation. Just run this code and you'll see:

a=[1 2 3 2 1];

b = zeros(1,11);

b(1:5) = a % Start at element 1

std1 = std(b)

b = zeros(1,11);

b(5:9) = a % start at element 5

std1 = std(b)

It's 1.078719779941187 for both cases. I invite you to give an example of where the placement affects the std dev.

And if the weighting can't be 0 then let it be "eps." That will minimize the std dev.

[Charles]

The original post had something missing. Here is an example

a =[1 2 3 2 1];

b = zeros(1,11);

b(1:5) = b(1:5) + a;

b(4:8) = b(4:8) + a * 1.5;

b(6:10) = b(6:10) + a * 2;

b(7 : 11) = b(7 : 11) + a * 2;

std(b)

c = zeros(1,11);

c(1:5) = c(1:5) + a;

c(3:7) = c(4:8) + a * 1.3;

c(6:10) = c(6:10) + a ;

c(7 : 11) = c(7 : 11) + a * 3;

std(c)

Answer 1

[Teja Muppirala]

Your problem can be expressed as an overdetermined system:

Given A and b where

A = [ 1 2 3 2 1 0 0 0 0 0 0; 
      0 1 2 3 2 1 0 0 0 0 0;
      0 0 1 2 3 2 1 0 0 0 0;
      0 0 0 1 2 3 2 1 0 0 0;
      0 0 0 0 1 2 3 2 1 0 0;
      0 0 0 0 0 1 2 3 2 1 0;
      0 0 0 0 0 0 1 2 3 2 1]

b =  [1 1 1 1 1 1 1 1 1 1 1];

Find an x (with size 1 by 7) with a maximum of 4 nonzero entries such that

norm(x*A - b) is minimized.

This can be done quite simply using plain old matrix division.

Step 0. Generate A and b.

A = [1 2 3 2 1 0 0 0 0 0 0];
A = gallery('circul', A);
A = A(1:7,:)

b = ones(1,11);

Step 1. Generate all unique combinations of 4 columns from 6 possible columns.

P = perms(1:7);
P = unique(sort(P(:,1:4),2),'rows');

Step 2. For each of those possible column combinations, use the pseudoinverse to find the optimal x, and the associated error for each case.

for n = 1:size(P,1)
  A4 = A(P(n,:),:);
  x{n} =  b/A4;
  err(n) = norm(x{n}*A4 - b);
end

Step 3. Find the combination that yielded the minimum error

[val, loc] = min(err);
x{loc} %<--- optimal weights
P(loc,:) %<--- optimal positions
plot(x{loc} * A(P(loc,:),:))

[Ivan van der Kroon]

You have my vote by putting it so straight-forward. However, as I understand it the positions are not unique (same values are allowed). So P will be a very large matrix if the number of placements is larger. Is there not another way to constrain it to integer values?

Typo: A misses a row; A=A(1:7,:);

[Teja Muppirala]

>> Typo: A misses a row; A=A(1:7,:);

Ah. Fixed that. Thanks.

If the scope of the problem is too much larger, then certainly trying to solve for every single case like what is shown here will not be computationally feasible. In that case you'd have to set it up as some sort of optimization problem. The genetic algorithm solver can be set up to deal with integer constraints.

[Teja Muppirala]

The fact that the positions can be non-unique is handled automatically. I am choosing 4 unique positions, but if there were a case where only using 3 were indeed more optimal than 4, then the linear least squares solution would indicate that by giving me a zero in one of the weights.

[Charles]

Hi Ivan and Teja, thank you guys for your comments and suggestions. Both of you have votes from me for your input. I will study your examples and hopefully, can make progress from as your suggestions don't quite solve my problem.

The vectors 'a' and 'b' are actually much larger ('a' has more than 100 elements).

[Ivan van der Kroon]

How large is b and are there any restrictions on the weights, e.g. nonnegative?

[Charles]

b is about 600 elements and the weights are non-negative.

Answer 2

[Ivan van der Kroon]

Allow me to define:

N=length(b);
M=length(a);
Q=3;%number if placements
p0=ones(Q,1);%starting vector for placements
w0=ones(Q,1);%starting vector for weigths

I reduced the p's to a single interger value such that you only need to give the first placement. For example

n=0:M-1;
b(n+p(j)) = b(n+p(j)) + a;

Now note that the entries of p are integers and must be within [1,2,...,N-M+1] otherwise the indices are out of bounds.

Now the problem is that I don't know how to specify that p should be these values in a standard optimization problem. If you can find this out (maybe just with ‘round’ if it is only of length 3 or 4, because it won’t cost you that much), I would suggest to use the optimization toolbox with the vector to look for is x=[p;w];

[Charles]

Hi, Ivan, thanks for your comment and I understand your point. Effectively, what you did is reduce the dimension of p's to a scalar. That point is appreciated.

Regarding this:

"Now the problem is that I don't know how to specify that p should be these values in a standard optimization problem." ==> Isn't that like specifying an upper and lower bound for the values of p and having the optimization algorithm find the optimal values of p's and w's?

BTW, I do not know how to solve the problem with optimization toolbox. Once again, thanks for your comment.

[Ivan van der Kroon]

So you have the known inputs a,b and Q. You want to minimize std(b); The point is that you have 2*Q variables, namely p and q, since they are both of length Q. Furthermore p is limited to only N-M+1 possible values whereas w can be varied continuously.

Hence you want to use a function like fminsearch to minimize std(b) as a function of x=[p;w].

Btw, did you try to solve this without numerical force?

[Charles]

Btw, did you try to solve this without numerical force? ==> I wish I knew how to.

[Ivan van der Kroon]

So what is the interval for w? And is std(b) normalized for this or must it always have the same mean? For instance, for Q=2, w=[1,2] is practically the same as w=[10,20]; If w is allowed to be negative the problem will be practically trivial (the triviality was already mentioned a while ago and I missed your response on it).

Without numerical force: minimizing std(b) is the same as minimizing Sum((b-mean(b))^2). Further, b(1)=a(1)*(p==1), b(2)=a(2)*(p==1)+a(1)*(p==2), etc (p==1 must be read as the number of times it equals the value 1). With weights it is b(1)=a(1)*( w(1)*(p(1)==1)+w(2)*(p(2)==1)+...+w(Q)*(p(Q)==1) ). Anyhow, you'll will yield length(b) equations which quadratic sum has to be minimized for Q or 2Q variables. Can you get somewhere from here?

[Ivan van der Kroon]

Another solution to the triviality will be to constrain Sum(w)=W for some predefined constant W.

[Charles]

Hi, thanks for your post. 'w' can only positive. I will go thru your posts.