I am looking for a reference where it is proven that the harmonic mean
$$\bar{x}^h = \frac{n}{\sum_{i=1}^n \frac{1}{x_i}}$$
minimizes ( in $z$) the sum of squared relative errors
$$\sum_{i=1}^n \left( \frac{(x_i - z)^2}{x_i}\right).$$
You could point out that this is a weighted least squares regression with weights $1/x_i$.
To make the connection with the references, revert to a standard notation in which you seek to find $\beta$ that minimizes $$\sum \omega_i (y_i - \beta)^2.$$
This is a model with a single constant regressor $$X = \pmatrix{1\\1\\\vdots\\1}$$ and weights matrix $$W = \pmatrix{\omega_1 & 0 & \cdots & 0 \\ 0 & \omega_2 & \cdots & 0 \\\vdots & \vdots & \ddots & 0 \\ 0 & \cdots & 0 & \omega_n}.$$
I have renamed "$x_i$" as "$y_i$" (the "response") and the parameter to be estimated is $\beta$ instead of $z$. The weights are $\omega_i=1/x_i$. It is necessary that they all exceed $0$. The solution is
$$\hat\beta = (X^\prime W X)^{-1}X^\prime W y = \frac{\sum_i x_i\omega_i }{\sum_i \omega_i} = \frac{\sum_i x_i/x_i }{\sum_i 1/x_i} = \frac{n}{\sum 1/x_i},$$
QED.
The same analysis applies to any positive sets of weights, providing a generalization of the harmonic mean and a useful way to characterize it.
When, as in a controlled experiment, the $x_i$ are viewed as fixed (and not random), the machinery of weighted least squares provides confidence intervals and prediction intervals, etc. In other words, casting the problem into this setting automatically gives you a way to assess the precision of the harmonic mean.
Viewing the harmonic mean as the solution to a weighted problem provides insight into its nature and, especially, to its sensitivity to the data. It is now clear that the most important contributors are those with the smallest values of $x_i$--and their importance has been quantified by the weights matrix $W$.
Douglas C. Montgomery, Elizabeth A. Peck, and G. Geoffrey Vining, Introduction to Linear Regression Analysis. Fifth Edition. J. Wiley, 2012. Section 5.5.2.
Why do you need a reference? This is a simple calculus problem: For the problem as you have formulated it to make sense, we must assume that all $x_i > 0$. Then define the function $$ f(z) = \sum_{i=1}^n \frac{(x_i - z)^2}{x_i} $$ Then calculate the derivative with respect to $z$: $$ f'(z) = -2\cdot \sum_{i=1}^n (1-\frac{z}{x_i}) $$ then solving the equation $f'(z)=0$ gives the solution. Now, of course we must check that this is indeed a minimum, for that calculate the second derivative: $$ f''(z) = -2 \cdot \sum_{i=1}^n (0-\frac1{x_i}) = 2 \sum_{i=1}^n \frac1{x_i} > 0 $$ for the last inequality we used, finally, that all $x_i>0$. Without that assumption we could indeed risk that we had found a maximum!
As for a reference, maybe https://en.wikipedia.org/wiki/Fr%C3%A9chet_mean or https://en.wikipedia.org/wiki/Harmonic_mean or references therein.