How can I find a functional square root in Mathematica?

I was not able to get a half iterate for $cos(...)$ anything, and from reading around a bit it appears that half iterates of $cos$ might be impossible either due to convergence or the evenness of $cos$'s series expansion terms. However, I was able to get a half-iterate for a small part of the domain of $sin(4 pi x)$ through fixed point iteration on the series, though it becomes inaccurate quite quickly:

(* Try to find a half iterate of Sin[4 [Pi] x] *)
halfit[x_] = Nest[(Sin[4 [Pi]*Normal[InverseSeries[Series[#, {x, 0, 6}]]]] + #)/2 &, x, 8];

Plot[{halfit[halfit[x]], Sin[4 [Pi] x]}, {x, -[Pi]/2, [Pi]/2}, 
 PlotRange -> {-1, 1}, 
 PlotStyle -> {Directive[Thick, Red], Directive[Blue]}]

I was able to get an approximation of the half-sine by a different method using a Newton series, though this does not work for a higher frequency sine like $sin(4 pi x)$ and produces a very noisy function. The resulting $mathrm{hsin}(mathrm{hsin}(x))approxsin(x)$ is not too bad an approximation judging by the plot:

newtonfhalf[f_, x_, mmax_] := 
 Sum[Binomial[1/2, m] Sum[
    Binomial[m, k] (-1)^(m - k) Nest[f, x, k], {k, 0, m}], {m, 0, mmax}]

nth = Function[{x}, newtonfhalf[Sin[#] &, x, 40]];
nthh2 = nth[nth[x]];
Plot[{Sin[x], nthh2}, {x, -4, 4}, 
 PlotStyle -> {Directive[Thick, Blue], Directive[Red]}]

---

I've had some luck with a neural network approach to the problem. I've found it's possible to train a network in a non-standard way to find an approximate half-iterate. Assume a network $N$ of 1 input and 1 output node with arbitrary layers in between and that we're trying to find a half-iterate for the function $mathrm{target}(x)$:

1. Run $N$ forward on a random input $x_i$ and generate output $y_i$
2. Run $N$ forward again using $y_i$ as input, generating output $y_i'$
3. The loss is $(mathrm{target}(x_i) - y_i')^2$. Back-propagate and update $N$ and return to step 1.

The resulting network is hopefully trained such that $N(N(x)) approx mathrm{target}(x)$.

I wasn't sure how to approach this in Mathematica but this is my first time using PyTorch ever, so what follows may be a bit basic:

import torch
import torch.nn as nn
import torch.optim as optim

from math import pi, sin, cos
import random
import csv


def targetfn(x):
    return sin(x)


class Net(nn.Module):

    def __init__(self):
        super(Net, self).__init__()
        self.lin = nn.Linear(1, 20)
        self.lmid1 = nn.Tanh()
        self.lmid2 = nn.Linear(20, 20)
        self.lmid3 = nn.Tanh()
        self.lout = nn.Linear(20, 1)

    def forward(self, w):
        w = self.lin(w)
        w = self.lmid1(w)
        w = self.lmid2(w)
        w = self.lmid3(w)
        return self.lout(w)


def train():
    net = Net()
    print(net)

    optimizer = optim.SGD(net.parameters(), lr=0.01)
    criterion = nn.MSELoss()

    # init random
    net.zero_grad()
    outinit = net(torch.randn(1))
    outinit.backward(torch.randn(1))

    for i in range(100000):
        x = random.uniform(-2 * pi, 2 * pi)
        target = torch.tensor([targetfn(x)])
        y1 = net(torch.tensor([x]))
        net.zero_grad()
        optimizer.zero_grad()
        y2 = net(y1)
        loss = criterion(y2, target)
        loss.backward()
        optimizer.step()

    return net


def main():
    net = train()

    with open("hfn.csv", 'w', newline='') as csvfile:
        csvwriter = csv.writer(csvfile, delimiter=',')
        n = 2000
        xmin = -2 * pi
        xmax = 2 * pi
        step = (xmax - xmin) / n
        x = xmin
        for i in range(n):
            csvwriter.writerow([x, net(torch.tensor([x])).item()])
            x += step


if __name__ == '__main__':
    main()

... and plotting in Mathematica:

data = Import["hfn.csv"];
intp = Interpolation[data];
Plot[{Sin[t], intp[intp[t]]}, {t, -2 [Pi], 2 [Pi]}, 
 PlotRange -> {-1.3, 1.3}, 
 PlotStyle -> {Directive[Thick, Blue], Directive[Thin, Red]}, 
 PlotTheme -> "Scientific"]

This is looking good for $sin(x)$. What about $cos(x)$? I changed targetfn in the python code above and at least I got something that looked close to a cosine wave: