Parameterization regression of rotation angle
The second way, predicting $x=cos(\alpha)$ and $y=sin(\alpha)$ is totally okay. 

Yes, the norm of the predicted $(x, y)$ vector is not guaranteed to be near $1$. But it is not likely to blow up, especially if you use sigmoid activation functions (which are bounded by they nature) and/or regularize your model well. Why should your model predict a large value, if all the training samples were in $[-1, 1]$?

Another side is vector $(x,y)$ too close to $(0,0)$. This may sometimes happen, and could indeed result in predicting wrong angles. But it may be seen as a benefit of your model - you can consider norm of $(x,y)$ as a measure of confidence of your model. Indeed, a norm close to 0 means that your model is not sure where the right direction is. 

Here is a small example in Python which shows that it is better to predict sin and cos, that to predict the angle directly:

# predicting the angle (in radians)
import numpy as np
from sklearn.neural_network import MLPRegressor
from sklearn.model_selection import cross_val_predict
from sklearn.metrics import r2_score
# generate toy data
np.random.seed(1)
X = np.random.normal(size=(100, 2))
y = np.arctan2(np.dot(X, [1,2]), np.dot(X, [3,0.4]))
# simple prediction
model = MLPRegressor(random_state=42, activation='tanh', max_iter=10000)
y_simple_pred = cross_val_predict(model, X, y)
# transformed prediction
joint = cross_val_predict(model, X, np.column_stack([np.sin(y), np.cos(y)]))
y_trig_pred = np.arctan2(joint[:,0], joint[:,1])
# compare
def align(y_true, y_pred):
    """ Add or remove 2*pi to predicted angle to minimize difference from GT"""
    y_pred = y_pred.copy()
    y_pred[y_true-y_pred &gt;  np.pi] += np.pi*2
    y_pred[y_true-y_pred &lt; -np.pi] -= np.pi*2
    return y_pred
print(r2_score(y, align(y, y_simple_pred))) # R^2 about 0.57
print(r2_score(y, align(y, y_trig_pred)))   # R^2 about 0.99


You can go on and plot the predictions, to see that predictions of the sine-cosine model are nearly correct, although may need some further calibration:

import matplotlib.pyplot as plt
plt.figure(figsize=(12, 3))
plt.subplot(1,4,1)
plt.scatter(X[:,0], X[:,1], c=y)
plt.title('Data (y=color)'); plt.xlabel('x1'); plt.ylabel('x2')
plt.subplot(1,4,2)
plt.scatter(y_simple_pred, y)
plt.title('Direct model'); plt.xlabel('prediction'); plt.ylabel('actual')
plt.subplot(1,4,3)
plt.scatter(y_trig_pred, y)
plt.title('Sine-cosine model'); plt.xlabel('prediction'); plt.ylabel('actual')
plt.subplot(1,4,4)
plt.scatter(joint[:,0], joint[:,1], s=5)
plt.title('Predicted sin and cos'); plt.xlabel('cos'); plt.ylabel('sin')
plt.tight_layout();




Update. A navigation engineer noticed that such a model would be most accurate when the angle is close to $\frac{\pi N}{2}$. Indeed, near 0° and 180° the angle $\alpha$ is almost linear in $\cos(\alpha)$, and near 90° and 270° it is almost linear in $\sin(\alpha)$. Thus, it could be beneficial to add two more outputs, like $z=\sin(\alpha+\frac{\pi}{4})$ and $w=\cos(\alpha+\frac{\pi}{4})$, to make model almost-linear near 45° and 135° respectively. In this case, however, restoring the original angle is not so obvious. 

The best solution may be to extract coordinates $(x,y)$ from both representations (in the second one, we need to rotate $(z,w)$ to get $(x,y)$), average them, and only then calculate arctan2.