Your summand is symmetric with respect to $k$ and $i$:

$$f(n,s) = \frac{1}{n}\sum_{k = 1}^n \sum_{i = 1}^{k} \bigg(\frac{\gcd(k,i)}{\operatorname{lcm}(k,i)}\bigg)^s$$

We can sum along skew diagonals to evaluate the sum. That is, we can convert $(k,i)$ to polar form $(\sqrt{k^2 + i^2}, tan^{-1}\frac{i}{k})$. By symmetry of the summand, the rays from the origin have the same value.

That is, when $\theta = tan^{-1}\frac{i}{k}$, $i = k \tan\theta$. We can vary $\theta$ from $[0, \frac{\pi}{4}]$. The $\gcd(j,j\tan \theta)$ is independent of $j$ but depends on $\theta$, we can therefore use $n$:

$$\frac{1}{n}\int_{0}^{\frac{\pi}{4}} \sum_{j=0}^{n} \left[\frac{\gcd(j,j\tan\theta)}{\operatorname{lcm}(j,j\tan\theta)}\right]^{s} d\theta = \int_{0}^{\frac{\pi}{4}} \left[\frac{\gcd(n,n\tan\theta)}{\operatorname{lcm}(n,n\tan\theta)}\right]^{s} d\theta = \sum_{k=0}^{n} \frac{1}{k^{s+1}} \rightarrow \zeta(s+1)$$

The polar integral(The integral is not continuous on $\theta$ but over the irrational($\tan\theta = \frac{i}{k}$ implies $\theta$ is irrational) values that correspond to the rays) to the sum is calculated easily enough because the values along each ray is constant w.r.t to $j$. Every ray corresponds to one of the values of $\frac{1}{k^s}$ but they have a weight of $\frac{1}{k}$.

E.g., the ray that accumulates $1$ has $\theta = \tan^{-1}(1) = \frac{\pi}{4}$, for $\frac{1}{2^s}$ it is $\theta = \tan^{-1}(2)$ but it has $\frac{1}{2}$ the density of $1$. Similarly for all the others.

If you are having trouble following this, simply look at $\frac{\operatorname{lcm}(k,i)}{\gcd(k,i)}$ in "polar" form, I'll make it easy for you(the text format obscures the patterns but they are there):

\begin{matrix}
\color{green}1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 & \\
\color{blue}2 & \color{green}1 & 6 & 2 & 10 & 3 & 14 & 4 & 18 & 5 & \\
\color{red}3 & 6 & \color{green}1 & 12 & 15 & 2 & 21 & 24 & 3 & 30 & \\
4 & \color{blue}2 & 12 & \color{green}1 & 20 & 6 & 28 & 2 & 36 & 10 & \\
5 & 10 & 15 & 20 & \color{green}1 & 30 & 35 & 40 & 45 & 2 & \\
6 & \color{red}3 & \color{blue}2 & 6 & 30 & \color{green}1 & 42 & 12 & 6 & 15 & \\
7 & 14 & 21 & 28 & 35 & 42 & \color{green}1 & 56 & 63 & 70 & \\
8 & 4 & 24 & \color{blue}2 & 40 & 12 & 56 & \color{green}1 & 72 & 20 & \\
9 & 18 & \color{red}3 & 36 & 45 & 6 & 63 & 72 & \color{green}1 & 90 & \\
10 & 5 & 30 & 10 & \color{blue}2 & 15 & 70 & 20 & 90 & \color{green}1 & \\
\end{matrix}

If you look you can see $k$th ray which has the constant value $\frac{\color{green}1}{k^s}$(displayed as just $k$ in the table) but they repeat at a rate of $\frac{\color{green}1}{k}$ along the ray.

Alternatively, if write the table in polar coordinates(we are rotating coordinate space 45 degree's) we get
\begin{matrix}
\color{green}1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \\
\color{green}1 & \color{blue}2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \\
\color{green}1 & 0 & \color{red}3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \\
\color{green}1 & \color{blue}2 & 0 & 4 & 0 & 0 & 0 & 0 & 0 & 0 & \\
\color{green}1 & 0 & 0 & 0 & 5 & 0 & 0 & 0 & 0 & 0 & \\
\color{green}1 & \color{blue}2 & \color{red}3 & 0 & 0 & 6 & 0 & 0 & 0 & 0 & \\
\color{green}1 & 0 & 0 & 0 & 0 & 0 & 7 & 0 & 0 & 0 & \\
\color{green}1 & \color{blue}2 & 0 & 4 & 0 & 0 & 0 & 8 & 0 & 0 & \\
\color{green}1 & 0 & \color{red}3 & 0 & 0 & 0 & 0 & 0 & 9 & 0 & \\
\color{green}1 & \color{blue}2 & 0 & 0 & 5 & 0 & 0 & 0 & 0 & \color{green}10 & \\
\end{matrix}

Where now the $k$th column is the "$k$th" ray. I.e., the first column in the above table corresponds to the diagonals/rays in the table above it.

actuallyclaim for the two sides to be equal, but rather them getting asymptotically closer? By your numerical data, a reasonable conjecture might be that $LHS-1\sim\zeta(s+1)-1$ in the sense of the ratio of the two tending to $1$. Since $\zeta(s+1)-1\sim 2^{-s-1}$, it should be enough to look at which terms in your LHS are $1$ and $2^{-s}$ and crudely estimating the rest. $\endgroup$ – Wojowu Dec 18 '19 at 10:57