compute a gradient with "minimal" colorstops from a ordered set of RGB colors I got a array of colors that are closely related = most of the time there is no sharp transition between color[n] and color[n+1].

Exemple of such an array of colors 
[0,.1,.2,.3,.4,.5,.6,.7,.8,.9,1.0,.8,.6,.4,0.2,0]

You'll notice that there are 3 "special" values. The colors
slowly increase from value 0 to 1 (in 10 steps) then
decrease more rapidly from 1 to 0 (in 5 steps).

So our array could also be described this way :  

an array of 16 colors where
Color[0]  whose value is  0 evolves linearly until
color[10] whose value is  1 evolves linearly to
color[15] whose value is  0

I don't really care about the length of my initial array. It will be  stretched (interpolated) so instead of using indexes it makes more sense to use percentage values
index 0  /15 = 0%
index 10 /15 = 66%
index 15 /15 = 100%  

so the initial array could be rewritten like this
[{0, 0%}, {1, 66%}, {0, 100%} ]

This last notation is close to what is called a gradient in svg or css. The special color are named "colorstops" and gradients are written this way
linear-gradient(black 0%, white 66%, black 100%);

One could notice that some gradients are equivalent  

linear-gradient(black 0%, white 100%) =
linear-gradient(black 0%, grey 50%, white 100%) =
linear-gradient(black 0%, lightgrey 25%, grey 50%, darkgrey 75%,  white 100%)

like 

interpolation[0, 1] =
interpolation[0, .5, 1] =
interpolation[0, .25, .5, .75, 1]

So my challenge is to go from
[0,.1,.2,.3,.4,.5,.6,.7,.8,.9,1,.8,.6,.4,.2,0]
to
linear-gradient(black 0%, white 66%, black 100%)
and not the naive
linear-gradient(
 black    0%,
 grey++++ 6%,
 grey+++  13%,
 grey++   20%,
 grey+    26%,
 grey     33%,
 grey-    40%,
 grey--   46%,
 grey---  53%,
 grey---- 60%,
 white    66%,
 grey---  83%,
 grey-    80%,
 grey+    86%,
 grey+++  92%,
 black   100%)  

The expected size of color array is no more that 2000 RGB colors and the goal for color-stops is no more that 10 and ideally 2.

EDIT
 in fact what i want to do is compute the blue lines like on the following graph
What you want to do is fit a piecewise linear function to your data.

No noise

If the data is perfect and has no noise, this is easy:


Form a line from the first points.  Remember those two points as the endpoints of the line.
If the third point lies on this line, discard the right endpoint and replace it with the third point.  Otherwise, start a new line from the second point and third point.
Repeat.


That can be done with a simple linear scan.

Noisy data

If the data is imperfect and noisy, you want to do a best fit of a piecewise linear function to your data.

There's no universal answer, as there are multiple ways to fit a piecewise linear function.  If you have $n$ points, you can always get a perfect fit by drawing $n-1$ lines (one line between each consecutive pair of points), but the resulting function might be very wiggly and not what you are looking for.  In practice, we probably have the idea that a solution with few pieces (few lines) is probably more likely than one with more pieces, all being equal.  On the other hand, with too few pieces, you might get a very poor fit.

Thus, you need a parameter to balance the number of pieces against the amount of error in your fit, and then you can try to find an optimal solution.  So, the first step will be to define a loss function that quantifies how good/bad each candidate solution is.  Once you have that, there are various algorithms for finding the best fit.  One approach is to use dynamic programming.  See https://en.wikipedia.org/wiki/Segmented_regression, https://dsp.stackexchange.com/q/1227/5874, https://stackoverflow.com/q/22028529/781723, and probably many other sources.

My guess is that if there is any noise/error in  your data, it will be very small, so any reasonable procedure will work well.