How to ensure anti-aliasing in image resampling at lower resolution?
A "1D image" is just a discrete signal. All the classic signal theory applies:


  Will such transformation introduce aliasing (and if so, what is the correct way to resample with averaging over pixels?)


So, what happens when you downsample is decimation. The aliases are the "folding" of existing spectrum over the conserved part of the spectrum; imagine you have a signal whose spectrum (for simplicity, assume that's its discrete Fourier transform) is 

$[A\,|\,B\,|\,C\,|\,D\,|\,E\,|\,F\,|\,G\,|\,H]$

and you simply throw away half of the samples, the resulting spectrum is

$\frac 12 [A+E| B+F| C+G| D+H]$. 


  Is there some math procedure, using which I can answer questions like #1 myself?


So what you'd need to do is suppress $[E\,|\,F\,|\,G\,|\,H]$ before throwing away the samples. Something that does that is usually called a low pass filter (LPF).

LPFs are often implemented as Finite Impulse Response Filter (FIR), where you take the current input sample, multiply it with a factor, add a previous input sample (also multiplied with a, possibly different, factor).

As it turns out, your averaging is indeed such an LPF FIR. It has the taps (coefficients) of 

$[0.5, 0.5]$. To see this, consider your $2i+1$th sample to be your "current" one, and the $2i$th to be "one in the past".

Now, such a FIR can be analyzed; it's z-domain frequency response is this:



Notice how at the vertical bar, i.e. the point below which you want to preserve and above you want to suppress the spectrum, this only has a value of 0.5; in fact, about 22% of the total pixel "energy" that passes through your filter is in the "upper" spectral half. That's too much, and you'll get aliasing.

You simply might want to build a better filter. This implies you'll add more taps, i.e. more samples right and left of your current sample

There's plenty of tools out there to design filters. Scipy's filter design tools will absolutely do; just design a low pass filter. 

For example, I've used GNU Radio's gr_filter_design tool, and specified the passband to end just below the half-band boundary (at frequency $0.45f_{sample}$), and the stopband to start just above that boundary ($0.55f_{sample}$)¹. Suppression should be 30dB, i.e. the magnitude in the $\le\frac1{1000}$ of the passband in the stopband. The resulting filter:



It has 13 taps, i.e. one coefficient for the "center" pixel and one for each of the six pixels left and right of it:



-0.00405147997662425, 0.009050002321600914, -0.02354922890663147, 0.04652633145451546, -0.07230162620544434, 0.09264583140611649, 0.9033603072166443, 0.09264583140611649, -0.07230162620544434, 0.04652633145451546, -0.02354922890663147, 0.009050002321600914, -0.00405147997662425


Notice the symmetry (typical for the class of linear-phase FIRs).

Original | Filtered with your averaging filter, then downsampled | filtered with the 13-tap filter, then downsampled (please view in full size; it doesn't make sense to compare things scaled by your browser)².





Visual comparison shows that the readability of the text suffers a bit more from the long filter response, while the reproduction of the fine structures of the dark cacao powder atop the tasty, juicy FFTiramisu is is better with the 13tap filter.
This is most probably an effect of the filter design (it's a design processed based on the Hamming window). Pre-scaling low pass filters employed by picture processing software often aren't even linear; that makes modeling significantly more complicated, but since sensory input in humans typically has logarithmic charactistics, I presume that hides "border effects" better.


  How, in general, discrete signals are anti-aliased? From DSP books I know that for continuous signals you just apply low-pass filter to your input signal and then sample whatever you got after the low-pass transform. But here is the problem: for images from the very beginning I have only discrete signal and I cannot figure out how to use continuous coordinate low-pass filters for discrete signals. May be I should treat an image as a piece-wise continuous function, which is constant within each interval [i,i+1)[i,i+1) ?


Well, in DSP downsampling is a very common thing to do, so you do with digital filters, as explained above. Don't consider your signal to be continuous; it's not.

¹ FIRs can't have infinitely steep responses. As a general rule, the closer the passband edge is to the stopband edge, the more taps you will need to realize the filter.

² Processing was done using Python, numpy, scipy.signal.lfilter, scipy.misc.im{read,save}.