What transformation can be used to establish relationship between two variables whose scatter plots looks like the plot below?
I have applied log Y, log X, exp Y, exp . But a relationship could not be established with high R-squared.
Note that the spread of $y$ increases linearly as $x$ increases.
You also seem to have (at least approximately) a linear increase $y$ as $x$ increases.
That is $E(Y|x) = \beta_0+\beta_1x$ but also $\sigma(Y|x) = \gamma_0+\gamma_1x$.
The conditional distribution doesn't seem to be particularly far from uniform (this might be important if you're interested in more efficient estimates, or if you're interested in prediction intervals).
There are approaches for the general case, such as iteratively reweighted least squares approaches; alternatively one might write down a suitable model for the conditional mean, variance and error distribution and solve by maximum likelihood, but some special cases (which seem to occur fairly often) are worth discussing.
If it were the case (impossible to tell from your plot) that the spread was proportional to $x$, i.e. that $\gamma_0$ was $0$, then we could proceed as follows:
That is $Y_i = \beta_0+\beta_1x_i+x_i\epsilon_i$, for $i=1,2,...,n$, where $\text{Var}(\epsilon)=\sigma^2I$.
Now divide through by $x_i$:
$Y_i/x_i = \beta_0 (1/x_i)+\beta_1+\epsilon_i$
and let $Y^*_i=Y_i/x_i$ and $x^*_i = 1/x_i$ we have:
$Y^*_i = \beta_1+\beta_0 x^*_i+\epsilon_i$
a constant variance-regression with the roles of the two coefficients swapped. When it applies, this is highly convenient.
Another possibility for this $\gamma_0=0$ case if all $y$ values are positive is to use a GLM with variance proportional to mean-squared (the gamma), or a similar constant coefficient of variance model (such as a Weibull regression). Again, the conditional distributions are more nearly uniform, so these regressions won't be as efficient as a choice closer to that. A third possibility when the variables are both positive and the regression intercept is also $0$ is to take logs of both sides; the relationship should remain linear and the conditional variance should become constant, but the error distribution might be quite skew. If the intercept on the regression line were not zero but was known, a shift could be applied before taking logs.
It won't matter much what method you use to deal with this; you will not get a considerably better fit than just sticking a like through the middle of the data (roughly, it looks like the intercept is just below 3 and the slope just above 1):
Any line that approximates $E(Y|x)$ for these data simply won't have a high $R^2$ (not that I'd pay much attention to $R^2$ in the best of circumstances, but even less so with heterosckedastic data).
Based on that $x$ feature alone, the relationship seems pretty clear: increasing conditional mean of $y$ as $x$ increases, with the conditional variance of $y$ also increasing as $x$ increases.
I do not see any transformations or expansions of the $x$ feature (e.g., polynomials or splines) that will do much for you.
You say that you want a high $R^2$, but I do not believe you to have the features to achieve that goal with any kind of generalizability (e.g., out of sample performance, adjusted $R^2$). You are missing some component that drives $y$.