Recommendations for Non-Proportional Hazards This is an issue that has plagued me for a long time and I have found no good answers in textbooks, Google, or Stack Exchange. 
I have data set of >100,000 patients for which four treatments are being compared. The research question is whether survival is different between these treatments after adjusting for a bunch of clinical/demographic variables. The unadjusted KM curve is below.

Non-proportional hazards were indicated by every method I used (e.g., unadjusted log-log survival curves as well as interactions with time and the correlation of Schoenfield residuals and ranked survival time, which were based on adjusted Cox PH models). The log-log survival curve is below. As you can see, the form of non-proportionality is a mess. Although none of the two-group comparisons would be too difficult to handle in isolation, the fact that I have six comparisons is really puzzling me. My guess is that I won't be able to handle everything in one model. 

I'm looking for recommendations on what to do with these data. Modeling these effects using an extended Cox model is likely impossible given the number of comparisons and differing forms of non-proportionality. Given that they are interested in treatment differences, an overall stratified model isn't an option because it won't allow me to estimate these differences. 
So, feel free to rip me apart, but I was thinking about initially estimating a stratified model to get the effects of the other covariates (testing the no-interaction assumption, of course), and then re-estimating separate multivariable Cox models for each two-group comparison (so, 6 total models). This way, I can address the form of non-proportionality for each two-group comparison and get a less wrong estimated HRs. I understand that the standard errors would be biased, but given the sample size, everything will likely be "statistically" significant.
 A: You certainly don't have marginal proportional hazards. That does not mean you don't have conditional proportional hazards! 
To explain in more depth, consider the following situation: let's suppose we have group 1, which is very homogeneous and has constant hazard = 1. Now in group two, we have a heterogeneous population; 50% are at lower risk than group 1 (hazard = 0.5) and the rest are at higher risk than group 1 (hazard = 3). Clearly, if we knew whether everyone in group 2 was a higher or lower risk subject, then everyone would have proportional hazards. This is the conditional hazards. 
But let's suppose we don't know (or ignore) whether someone in group 2 is at high or low risk. Then the marginal distribution for them is that of a mixture model: 50% chance they have hazard = 0.5, 50% they have hazard = 3. Below, I provide some R-code along with a plot of the two hazards. 
# Function for computing the hazards from 
# a 50/50 heterogenious population
mix_hazard <- function(x, hzd1 = 0.5, hzd2 = 3){
  x_dens <- 0.5 * dexp(x, hzd1) + 0.5 * dexp(x, hzd2)
  x_s    <- 1 - ( 0.5 * pexp(x, hzd1) + 0.5 * pexp(x, hzd2)) 
  hzd    <- x_dens/x_s
  return(hzd)
}

x <- 0:100/20
plot(x, mix_hazard(x), 
     type = 'l',
     col = 'purple', ylim = c(0, 2), 
     xlab = 'Time', 
     ylab = 'Hazard', 
     lwd = 2)
lines(x, rep(1, length(x)), col = 'red', lwd = 2)

legend('topright', 
       legend = c('Homogeneous',
                  'Heterogeneous'), 
       lwd = 2,
       col = c('red', 'purple'))


We see clearly non-proportional marginal hazards! But note that if we knew whether the subjects in group 2 were high risk or low risk subjects, we would have proportional hazards. 
So how does this affect you? Well, you mentioned you have a lot of other covariates about these subjects. It is very possible that when we ignore these covariates, the hazards are non-proportional, but after adjusting for them, you may capture the causes of the heterogeneity in the different groups, and fix up your non-proportional hazards issue. 
A: Fantastic question fantastic answers.  I'll add that you should consider a model making much different assumptions such as the lognormal survival model. Use the normal inverse function for the y_axis instead of log-log.   Still need to covariate adjust.  So also look at normality of residuals stratified by treatment. This is covered in a case study near the end of my course notes at https://hbiostat.org/rmsc
