Algorithm for ordering a list under a similarity function?
The variant of your problem where you DON'T square the distance is called the (weighted) Minimum Linear Arrangement problem. It's NP-hard and there are a number of different approximation methods for it. Here's one link that popped up that has code that you might find useful. 

But fortunately for you, you squared the distance, which changes things quite a bit. Let's expand the expression: 
$$
\begin{align*} 
\sum_{a,b \in A} s(a,b) (i(a) - i(b))^2 &amp;= \sum_{a,b \in A} s(a,b) (i^2(a) +i^2(b) - 2 i(a) i(b)) \\
&amp;= \sum_a i^2(a) \sum_b s(a,b) + \sum_b i^2(b) s(a,b) - 2 \sum_{a,b} s(a,b)i(a)i(b) \\
&amp;= 2 \sum_a S(a) i^2(a) - 2 \sum_{a,b} s(a,b)i(a)i(b) \\
&amp;= 2 u^\top A u 
\end{align*}
$$

where $u$ is the vector $i(a)$ and $A_{a,b} = -s(a,b), A(a,a) = \sum_b s(a,b)$.

All of this is a long-winded way of saying that you've essentially defined the Laplacian operator for the complete graph with weights given by your similarity function. At this point, minimizing the form is straightforward. You compute the decomposition of A in the form $A = V \Lambda V^T$, at which point your goal is to minimize the norm $\| u^T V \sqrt{\Lambda}\|$. But note that since $A$ is symmetric, $V$ is orthonormal, and is essentially a rotation matrix, and so what you really want to do is line up $u$ in reverse order of the eigenvalues of $A$ (the elements of $\Lambda$)