Is there a way to decompose a vector into orthogonal vectors using regression? I have a set of data Y. Each Y is a vector:
Yi = [Yi,1, Yi,2,...Yi,n].
I know that each Y is a mixture of A and B. I also know that A and B are orthogonal vectors, meaning that A . B = 0. For each i,  Yi = αi A + βi B + εi. Where αi, βi, B and A are all unknown. 
Are there any existing method that allow me to find those unknown coefficients (αi, βi) and vectors (A , B) that can minimize Σ(εi)2?
 A: Restatement of the problem
Consider each $Y_i$ to be a column vector with $n$ components and let $Y$ be the matrix whose columns are $Y_1, Y_2, \ldots, Y_k$ in any order.  Let $U$ be an $n\times p$ matrix: We have in mind $p=2$ and are thinking of the columns of $U$, say $(U_1, U_2, \ldots, U_p)$, as a basis for a subspace to which all the $Y_i$ are "close."  This would mean there exist $p$-vectors $\omega^{(i)}=(\omega^{(i)}_1, \omega^{(i)}_2, \ldots, \omega^{(i)}_p)^\prime$ for which the differences
$$\varepsilon^{(i)} = Y_i - U\omega^{(i)} = Y_i - (\omega^{(i)}_1 U_1 + \cdots + \omega^{(i)}_p U_p)$$
tend to be "small;" specifically, their sum of squares should be minimized.  If we assemble the $\omega^{(i)}$ into the columns of a $p\times k$ matrix $W$, we can express this criterion as minimizing the value of 
$$||Y - UW||_F^2$$
where the squared Frobenius norm $$||A||_F^2$$ of any matrix $A$ is the sum of squares of its components.  Since the rank of $U$ obviously does not exceed the number of its columns $p$, $UW$ is a minimum-norm rank-$p$ approximation of $Y$.
Analysis
The Frobenius norm is unchanged by right- and left-multiplication by orthogonal matrices (practically by the definition of orthogonal).  The Singular Value Decomposition, or SVD, of any matrix like $Y$ consists of applying such multiplications to reduce the matrix to a diagonal matrix $\Sigma$ (necessarily of the same dimensions as $Y$).  ("Diagonal" means that the only nonzero entries of $\Sigma$ must have equal row and column indexes.) Moreover, since row and column permutations as well as negating entire rows or columns are also orthogonal multiplications, it can be arranged that (1) all the entries of $\Sigma$ are nonnegative and (2) the values along the diagonal are decreasing: $$\Sigma_{11} \ge \Sigma_{22} \ge \cdots \ge \Sigma_{\min(k,n),\min(k,n)} \ge 0.$$
It is not hard to show that the minimum-norm rank-$p$ approximation of $\Sigma$ is the matrix $\Sigma_p^{*}$ obtained by replacing all diagonal entries beyond row $p$ with zeros.
Solution
Let, then,
$$Y = U \Sigma V^\prime$$
be the SVD of $Y$ where the columns of $U$ are mutually orthogonal.  As we have seen, for any $p$ not exceeding the number of columns of $Y$ (namely $k$),
$$Y^{*} = U \Sigma_p^{*} V^\prime$$
must be the best rank-$p$ approximation to $Y$ (in the Frobenius norm) because $\Sigma_p^{*}$ is the best rank-$p$ approximation to $\Sigma$.  Notice that
$$W = \Sigma_p^{*} V^\prime$$
has at most $p$ non-zero rows.  The first $p$ columns of $U$ and the $k$ column vectors of $W$ (stopping at row $p$) solve the problem: the former play the role of the orthogonal regressors $A$ and $B$ in the question while the latter are the coefficients obtained from regressing the columns of $Y$ on $A$ and $B$.
Specific solution
With $p=2$, take $A$ to be the first column of $U$ and $B$ to be the second column.  By construction, $A$ and $B$ will be orthogonal (and each will have unit length).  The $\alpha_i$ and $\beta_i$ can be read off the columns of $W$ (the $\alpha_i$ are in its first row and the $\beta_i$ are in its second row).
Here is the result of performing this calculation 5000 times for a set of random $12\times 6$ matrices of rank $2$, to which iid Normal errors with standard deviation $0.1$ were added.  This standard deviation can be estimated by dividing the squared Frobenius norm of the residuals by the degrees of freedom and taking the square root.  The degrees of freedom will be the amount of data in $Y$, $n\times k$, minus the number of parameters estimated.  The number of parameters needed to specify $p$ orthogonal vectors is $n + (n-1) + \cdots + (n-p+1) = (n+1)n/2 - (n+1-p)(n-p)/2$.  Another $pk$ parameters are estimated in $W$.  However, $p$ of them are redundant because the $p$ vectors have all been normalized to unit length.  Thus, we may compare
$$\frac{||Y - Y^{*}||_F^2}{nk - ((n+1)n/2 - (n+1-p)(n-p)/2 + pk - p)}$$
to $\sigma^2$ as a test of how well this procedure works.  Their square roots are plotted in the histogram.  For reference, the vertical red line indicates the value of $\sigma$ itself.

The comparison between simulation and theory is excellent.  It will become poorer when the magnitude of $\sigma$ becomes large compared to the average size of the components of $Y$, because then this procedure does not estimate $\sigma$ very well: substantial overfitting can occur.  But that's not of concern for the solution of the original problem.

The following R code produced the figure.  It contains an efficient implementation of the solution in the four lines following the "Find its rank-p approximation" comment.
sigma <- 0.1 # SD of the errors
n <- 12 # Number of rows of `Y`
k <- 6 # Number of columns of `Y`
p <- 2 # Rank--must not exceed min(k, n)
set.seed(17)
Norm.squared <- replicate(5e3, {
  #
  # Generate a random matrix Y by adding iid error terms to `p` vectors.
  #
  U <- matrix(rexp(n*k), ncol=k)
  W <- matrix(runif(k*k), ncol=k)
  if (p < k) U[, (p+1):k] <- 0
  Y <- U %*% W + rnorm(n*k, sd=sigma)
  #
  # Find its rank-p approximation.
  #
  X <- svd(Y)
  U.hat <- X$u[, 1:p]
  W.hat <- diag(X$d[1:p], p, p) %*% t(X$v[, 1:p, drop=FALSE])
  Y.hat <- U.hat %*% W.hat
  #
  # Report the squared Frobenius norm of the residuals.
  # Suitably adjusted, it should be close to `sigma`, 
  # because it estimates `sigma`.
  #
  epsilon.hat <- Y - Y.hat
  sum(zapsmall(c(1, epsilon.hat))[-1]^2)
})
df <- n*k - ((n+1)*n/2 - (n+1-p)*(n-p)/2 - p + p*k)
hist(sqrt(Norm.squared/df), breaks=35, xlab="Value")
abline(v=sigma, lwd=3, lty=1, col="Red")

