How to determine the cut off value of an hyperellipsoid in order to retrieve a single quantile of a multivariate normal distribution? Introduction
My goal is to retrieve the $\alpha$ quantile of a N(0, H) (multivariate normal) random variable $X$ where H is a known d-dimensional positive definite matrix (with $d >3$). In other words, my goal is to retrieve a single Value-at-Risk (VaR) of X at the $\alpha \in (0,1)$  confidence level, i.e.
$$\text{VaR}_{\alpha}(X) = \text{inf} \ \{x \in \mathbb{R}: F_X(x) > \alpha\}$$
I aim to retrieve the $\text{VaR}_\alpha$ for $\alpha \in \{0.1, 0.05, 0.01\}$ in the context of a portfolio comprising $d>3$ financial log returns where we assume that the portfolio log return follows a multivariate normal distribution, i.e. N(0, w^T H w). We define w as an equally weighted vector of length $d$ meaning that each asset is equally represented in our portfolio. In this context VaR is the maximum loss with 1-$\alpha$ probability. Given the definition of VaR, VaR is expected to increase with a decrease in the $\alpha$. The idea is to use VaR estimates in a backtesting analysis in order to evaluate different models used for forecasting the covariance matrix H.      
Approach
In case of unimodal symmetric distributions such as the multivariate normal distribution considered in the question, quantiles correspond to 
 Highest Density Regions. As pointed out here, the highest density region of an N(0,H) random variable is an ellipsoid centered at its mean, 0, and oriented per the covariance matrix H: 
$$x: x^TH^{-1}x \le y$$  
In order to find a single quantile that satisfies the definition of Value-at-Risk we proceed as follows: the set of solutions that satisfy the condition of the ellipsoid equation, the level set, is retrieved by randomly generating many points on this ellipsoid. Then, I evaluate the portfolio return at each solution satisfying the condition: for example, let's suppose $s = (x_1, x_2, \dots, x_d)$ is a solution/ point on the ellipsoid, then the portfolio return at this point will be $w\cdot s$, which corresponds to the mean of the solution vector because we stated that each asset is equally represented in the portfolio and $x_i \in s$ corresponds to the return of asset $i$. Consequently, the lowest return, the minimum, on the ellipsoid is the maximum loss. 
Question
Now, in order to ensure that the maximum loss on the ellipsoid is the maximum loss with 1-$\alpha$ probability, hence is the Value-at-Risk of X at the $\alpha$ confidence level, I need to define the cutoff value for the ellipsoid, i.e. $y$. At this point, I am a bit confused concerning the definition of the cutoff value for the ellipsoid.
Here it is stated that the cutoff value for the ellipsoid can be determined from the Chi-square with d degrees of freedom. For instance, the highest density region capturing 0.95 probability of the N(0,H) is found with y= value such that Chi-square with d degrees of freedom 0.95. When defining the cutoff value from the chi-square with d degrees of freedom and following the procedure explained in the approach section, I would interpret the lowest return on the ellipsoid as the maximum loss with 0.95 probability, i.e. $\text{VaR}_{0.05}$. Is it correct to state that these maximum loss values are then the maximum loss with 1-α probability, hence is the Value-at-Risk? VaR does increase with a decrease in $\alpha$ but at the same time VaR seems a bit a large for the context which makes me doubt my interpretation and definition of the cutoff value for the ellipsoid determined from the Chi-square with d degrees of freedom and its link to the $\alpha$ quantile of the distribution (Value-at-Risk). 
Any comments that could point me in the right direction are more than welcome.               
Python code for obtaining random points on the ellipsoid
import numpy as np
from scipy.linalg import sqrtm
from scipy.stats import ortho_group
from scipy.stats._continuous_distns import chi2

dim = 10
# Create a positive definite matrix H with shape (dim, dim).
# evals is the vector of eigenvalues of H.  This can be replaced
# with any vector of length `dim` containing positive values.
evals = (np.arange(1, dim + 1)/2)**2
# Use a random orthogonal matrix to generate H.
R = ortho_group.rvs(dim)
H = R.T.dot(np.diag(evals).dot(R))

# y determines the level set to be computed.
y = chi2.ppf(q=0.95, df=dim)    
# Generate random points on the ellipsoid
nsample = 100000
r = np.random.randn(H.shape[0], nsample)
# u contains random points on the hypersphere with radius 1.
u = r / np.linalg.norm(r, axis=0)
xrandom = sqrtm(H).dot(np.sqrt(y)*u)
# Compute maximum loss on the ellipsoid
xrandom_min = np.min(np.array([np.mean(x) for x in xrandom.T]))
print("Maximum loss i.e. Value-at-Risk:")
print(xrandom_min)

 A: It seems that you made an error very early in your story when you say

where we assume that the portfolio log return follows a multivariate normal distribution, i.e. N(0, w^T H w).

That distribution is a univariate normal distribution and not a multivariate normal distribution.

What you should use is
If $R=w^T x$ is your loss/return, then this means that you are looking for the distribution of a variable that is the linear combination of the components of the multivariate normal distributed variable
The sum of any linear combination of the components of a multivariate distributed normal distribution is a univariate normal distributed variable (it is one way how you can define a multivariate normal distribution).
In your case:
$$R \sim N(0, w^THw)$$
for background about this $w^THw$ term see Matrix notation for the variance of a linear combination
Now you can just look for the percentiles of this univariate normal distribution.
So you look for $\text{VaR}_\alpha(w^TX)=\text{VaR}_\alpha(R)$ instead of $\text{VaR}_\alpha(X)$

Why and how you get a difference:

*

*You need to look at the minimum of the 95% of samples that optimize the loss.

*Instead you use the minimum of the 95% of samples that optimize the Mahalanobis distance.

These two will cover two entirely different regions, or at least shapes (your loss isosurfaces are defined by $w^T x=c$, which are hyperplanes).
Thus your Mahalanobis ellipsoids should contain some points outside the region of the 95% samples with optimized loss. Therefore the calculated value in this way is lower than the VaR. Because inside the 95% ellipsoids you include some of 5% point beyond the 5% hyperplane with lowest cost.
Graphical example (comparing two 95% boundaries)
Below is an image that demonstrates the principle (in two dimensions but the principle remains the same for more dimensions). Thousand datapoints are simulated and drawn for $H=\begin{bmatrix} 4 & 1 \\ 1 & 2 \end{bmatrix}$. The colours depict varying return which is calculated as $w \cdot x$, and with equal weights this is $x_1+x_2$.
You see that the the boundary around the 95% points with highest density, is something different than the boundary for the 95% with the highest return (or lowest loss).

A: If you're just after VaR, one way to calculate VaR is the "Delta Normal" method. 
If $\Phi^{-1}$ is the inverse CDF function for a standard normal random variable, $\mathbf{w}$ is your nonrandom weight vector, and $\mathbf{H}$ is the covariance matrix of the future returns, then 
$$
\text{VaR}_{\alpha}(X) = \Phi^{-1}(\alpha)\sqrt{\mathbf{w}^T \mathbf{H} \mathbf{w}}.
$$
You can see the probability of your return being less than this lower bound is $\alpha$:
\begin{align*}
P\left(R \le \Phi^{-1}(\alpha)\sqrt{\mathbf{w}^T\mathbf{H} \mathbf{w}}\right) &= P\left(R(\mathbf{w}^T\mathbf{H} \mathbf{w})^{-1/2} \le \Phi^{-1}(\alpha)\right)\\
&= \Phi \left( \Phi^{-1}(\alpha)\right) = \alpha.
\end{align*}
This does not require simulation, and the only place where matrices and vectors come into play is when you take the variance of the scalar return: $\operatorname{Var}(R) = \operatorname{Var}(\mathbf{w}^T \mathbf{X}) = \mathbf{w}^T \operatorname{Var}(\mathbf{X}) \mathbf{w}= \mathbf{w}^T \mathbf{H} \mathbf{w}$.
Here's some code:
import numpy as np
from scipy.stats import norm

def compute_var(alpha, cov_mat, weight_arr):
    variance = np.dot(np.dot(weight_arr, cov_mat),weight_arr)
    std_dev = np.sqrt(variance)
    return norm.ppf(alpha)*std_dev

weight = np.repeat(1,3)/3.0 #equal weights
Sigma = np.array([[1.0,.5, .5],[.5,1.0,.5],[.5,.5,1.0]])
compute_var(.05, Sigma, weight)

