How do I split a vector and minimize absolute error efficiently? I have a vector $y\in\mathbb{R}^n$. I want to find a computationally efficient way to find a split point $k$, such that:
$\sum_1^k |y_i-\theta_1|+\sum_{k+1}^n |y_i-\theta_2|$ is minimised. The $\theta$ values are the median of the set they are in.
I can easily find such a point iterating amongst all possible $k$ values, but I assume there is a faster way to do it and a theoretical explanation (in the end this is just a simple case of change point detection).
I post a code sample that performs the task using a for loop in R. I do not care about the programming language for the answer, I would prefer to understand how it is done faster.
set.seed(1)
n<-20L
y<-rnorm(n)
besterror<-Inf
for(i in 1:(n-1L)){
  y1<-y[1:i]
  y2<-y[((i+1L):n)]
  theta1<-median(y1)
  theta2<-median(y2)
  error<-sum(abs(y1-theta1))+sum(abs(y2-theta2))
  if(error<besterror){
    k<-i
    besterror<-error
  }

}
print(paste("Best split at",k,"absolute error",besterror))

 A: I think it can be done in $O(n)$ time if $y$ is sorted. In my pseudo-R it would look something like
P = prefix_sums(y) // O(n)
S = suffix_sums(y) // O(n)
min_idx = -1
min_sum = Inf
low_idx, high_idx = 0, 0
low_val, high_val = 0, 0
for (i in 1:length(y)) {
    low_idx, low_val = med_low(y, i) // get median index and value in lower part
                                     // O(1) since y - sorted
    high_idx, high_val = med_high(y, i) // O(1)
    sum_low = low_idx * low_val - P[low_idx] 
            + (P[i] - P[low_idx - 1]) - (i - low_idx) * low_val 
    sum_high = (high_idx - i - 1) * high_val - (S[i + 1] - S[high_idx])
             + S[high_idx] - (n - high_idx + 1) * high_val
    if (sum_low + sum_high < min_sum)
        // update ...
}

Basically the idea is that we can calculate the two sums in constant time if we precomputed all partial prefix and suffix sums of vector $y$. Also, finding median in sorted vector is also $O(1)$
A: The changepoint.np R package does this in O(n) time using a range of quantiles.  
library(changepoint.np)
set.seed(1)
n<-20L
y<-rnorm(n)
out<-cpt.np(y)
cpts(out)

This uses the defaults and identifies that there are no changes in the data.  If you want to know which location is the most likely for a changepoint then you can use the CROPS penalty:
out1<-cpt.np(y,penalty="CROPS",pen.value=c(3,7))
cpts.full(out1)

The CROPS penalty gives you all segmentations within a range of penalties, using cpts.full we can see all the segmentations with a penalty between 3 and 7 in the above.  It gives from no changepoints at a penalty of 7 to eight changes for a penalty of 3 and all in between.
You can set nquantiles=1 but i'm not sure what quantile that would calculate, you would need to read the original paper here (open access):
PAPER
It isn't quite the answer to the question as it isn't a single changepoint, but it is a more general solution, potentially multiple changes and different quantiles, with a faster O(n) computational time than those listed already.
