Vect tip of the week : Vectorization of piecewise functions

Hi all.

As every week*, we will have a tip for those who want to vectorize their matlab programs, for reasons already stated in earlier posts. In this case, we will look at the vectorization of piecewise functions.

Let’s say you have a function like this :

f(x)=\left\{\begin{array}{lr}x,&x\in[0,\frac{1}{2}]\\x^2,&x\in(\frac{1}{2},1]\end{array}\right.

which, when plotted, would look like this

This is the graph of the piecewise function

OK. So, what’s the deal? Well, from the earlier posts on vectorization, you may have seen that we always did simple stuff. If you have coded a piecewise function before in Matlab, you know that it’s something like this

function y = pieceWise(x)
if(x>a1 && xa2 && x
   y = a different value here
.
.
.
else
  y = a value in the case all else fails!
end

So, this code would implement a piecewise function which changes in the intervals [a_{1},b_{1}],[a_{2},b_{2}]\ldots. For our own example, the code would be

function y = piecewise2(x)
 if(0<=x && x<1/2)
 y = x;
 elseif(1/2<=x && x<1)
 y = x.^2;
 else
 y = 0;
 end
end

where we supposed that our function is zero outside of [0,1]. Our goal is to write that function in a format where if x is a vector then y = f(x), another vector.

Of course, trying that with the way our function is written will not work. Give that function a vector as input and you will get back

??? Operands to the ||
and && operators must
be convertible to
logical scalar values.
Error in ==> piecewise2
at 2
 if(0<=x && x<1/2)

That’s because the && operator is not built to handle vectors of logical expressions.

Wait, what? …  What’s a vector of logical expressions in Matlab?

Well, think of it this way. When we have a logical expression like x>0, then this expression can give back two possible numbers. Feed it with a positive x and it gives back 1.  Else it would just return 0.  But what would happen if we changed the scalar x with a vector \bar{x} of real numbers? Isn’t it logical to expect another vector with 0s and 1s, according to which elements of \bar{x} are greater (or smaller) than 0?

Either way, that’s what happens. If x=(0,10,-2,-3.1), then (x>0)=(0,1,0,0).

So, here’s a question for extra credit. Say we have y_{1}=(x>0)=(0,1,0) and y_{2}=(x<1/2)=(1,1,0). Then, what does the piecewise product of the two vectors, y_1 and y_2 stand for? I will wait for a little bit, so you can figure this out.

. . . . .

If you’ve thought about it, you found out that y_{1}.*y_{2}=(0,1,0) is nothing else than the vector of the conditional expression that would be generated by (x>0) \text{ AND } (x<1/2) (if that were possible).

So how do we use this idea? We can’t just go in our code and change it to

function y = piecewise2(x)
 if((0<=x).*(x<1/2))
 y = x;
 elseif((1/2<=x).*(x<1))
 y = x.^2;
 else
 y = 0;
 end
end

The problem here is that the if-clause is not gonna figure it out that you want to do a separate “if” for every element of the vector (x>0)\text{ AND }(x<1/2). It’s gonna only return the first value that it finds. Try it and see for yourselves.

So, how can we use this information? We only need to think this function from a different perspective.

Isn’t this function

f(x)=\left\{\begin{array}{lr}x,&x\in[0,\frac{1}{2}]\\x^2,&x\in(\frac{1}{2},1]\end{array}\right.

equal to this

f(x)=x\cdot X_{[0,1/2]}(x)+x^2\cdot X_{(1/2,1]}(x) ???

And isn’t the conditional expression x\geq 0 \text{ AND } x\leq 1/2 just X_{[0,1/2]}(x)? Note that X_{A}(x) is called the characteristic function of set A and it equals 1 if x belongs in that set, 0 if not.

Using all of those ideas, we can write our function as

function y = piecewise2(x)
y = ((0<=x).*(x<1/2)).*x + ((1/2<=x).*(x<1)).*x.^2;

Congrats! This function is now vectorized! With the same idea, you can vectorize complex conditionals statements by just thinking about the arithmetic between them. Assign values of 1 and 0 to the true and false statements and think how you could gain the same results.

Try to figure out how (x<-1/2)\text{ OR }(x>1/2) would work by doing arithmetic between those vectors.

Until next time.

*Well,  I am not very consistent with that “every week” promise. I should probably rename these posts as “Vect tip of the month”. 😛

Advertisements