Rather than an answer, I have a few questions which are relevant to finding efficiencies  for solving your problems.  You say that one of the complications is that you have many such systems.  So my first question is:  Does the coefficient matrix A vary from problem to problem?  Secondly, is the known matrix A of full column rank?  And finally, is there any additional information you have about the data  b such as error estimates or distribution?  This would be useful information in designing an algorithm.

Are you seeking the least-squares relative error (1+/-0.1, 100+/-10) rather than absolute error (1+/-10, 100+/-10), which would result from minimizing the residual in log(b) rather than b? That will work. If the only problem is the zero (or negative) values, then the log of (b+epsilon) or abs(b+epsilon) will work. You then adjust epsilon until you get a result that doesn't over emphasize the zero values. Also, you don't have to use log(b). You can use some other function, such as b/(a+b) or tanh(a*b). You also don't need to resort to a QR decomposition. The coefficients are given by multiply[ inverse[ multiply[ transpose[A],A]], multiply[ transpose[A],B]], where [A] is the equations and [B] are the solutions. You can easily implement this in an Excel spreadsheet, plot the results and watch how the comparison changes with epsilon and/or a. You can use the matrix formulas or the LINEST function.

Thank you for the answer, Dudley.

1) In general, I seek the "absolute" least-squares solution, i.e. minimization of the norm of the residual. Reformulating the problem in the relative sense may help in some situations, but I will likely have the problems with zero denominator.

2) (b+epsilon) is not a good way out because "not-overemphasizing the zero values" is not a quantifiable thing. And I cannot afford manual reviewing of the results, I need an automatic method.

3) The logarithm link function is motivated by the physics of the problem, and so it should not be changed to a rational function or tanh.

4) QR decomposition, in my mind, is the easiest way to solve an overdetermined linear system in the least-squares sense (because it naturally splits the system into a solvable triangle system and unsolvable "perpendicular" part). Most stock functions for solving overdetermined systems are based on QR decomposition.

After some thinking, I have realized that my problem is equivalent to multinomial logistic regression. And the Wikipedia's article on it states that the problem is usually solved (i.e. fitted) by iterative quasi-newtonian methods.

So it seems like there is no semi-analytical way to solve the problem fater.

You can accomplish this in an Excel spreadsheet using the Solver function. Begin with some estimate of the coefficients for the curve fit and then minimize some calculated residual. The Solver function uses a constrained minimization algorithm, so you can also limit the result to positive values.

I would pose this as a generalized linear model where you are not required to transform the data, but rather posit the model as log(E(b))=Ax. The solution is obtained by maximizing the likelihood under the assumption that the data "b" are from a lognormal, gamma, poisson (for count data), or other similar right skewed distribution.  

Note it is important to recognize that the left hand side is notational and does not imply any transformation of the actual data.  This model simply assumes that the log-mean is a linear function of the covariates--in this situation, least squares is not really appropriate for optimizing the model paramters due to heteroscedasticity.  

See McCullough and Nelder for an introduction and use R--function glm, to fit the models.

http://www.amazon.com/Generalized-Edition-Monographs-Statistics-Probability/dp/0412317605

Just noticed you noting that the data are binary...so yes you want logistic regression which is also a generalized linear model but with logit link function log(p/(1-p)) = Ax.  Other link functions are also possible for binary data such as complementary log-log. 

John, log(b)=A*x transformation is not possible because I have zero elements in b. And solving log(b)=A*x as a linear system (i.e., in the "log(b)-norm") results in oversensitivity in some cases: for example, if b=[1000,500,0] and the model predicts b*=[1000,500,1] then in the "log-norm" the norm of the residual is r=norm(log(b)-log(b*))=[+infinity] and, thus, such difference is extremely penalized, although from both common sense viewpoint and statistical viewpoint these vectors b and b* are "very close".

To be more precise, what I am doing is fitting multi-source capture-recapture models to data. From the statisitical point of view, I have multinomially distributed data, with one of bin counts being unmeasurable (and the rest measurable bin counts constitute vector b). The bin counts (effectively, the product of the size of the sample (unknown!) and multinomial probabilities) are predicted by term exp(A*x), where matrix A is defined by the "physics" of the model and x is the vector of hidden explanatory variables (mostly capture coverage probabilities and interactions).

I do understand that the most correct way to fit such a model is to use maximum likelihood approach. But to make the process of numerical optimization faster, I want to have a good initial approximation to the model parameters, and I plan to extract it from the (hopefully, fast) least-squares fitting of the model.

If the bottleneck is getting a good initial approximation for the nonlinear problem, one approach is to use the solution from a linearized variant that first deals with the zeros. To this end you might proceed as follows.

(1) Drop zeros from b. Call resulting data vector b'.

(2) Drop corresponding rows of A. Call resulting matrix A'.

(3) Find the usual least-squares solution to log(b') = A'x. Call solution x_0.

(4) Use x_0 as your initial approximation in the original setup, for a nonlinear least squares fit.

If b has "many" zeros this is not likely to give a good starting value, and in such cases you might want to consider using an epsilon, as was already suggested. This is only to get an initial approximation, and it might turn out that the problem is not terribly sensitive to the choice of epsilon for this particular purpose (that of course will require experimentation on your part to determine).

Daniel, I have already tried dropping zeros together with corresponding rows from A: the problem appears to be very sensitive to it. So for now I stick to effectively the "epsilon-approach": I replace zeros with ones to get an initial approximation (ones are effectively epsilons because most other counts in b are dozens and hundreds).

I see what you are needing, but just one point....the glm approach does not require log transformation of the dependent variable.  That is why I suggested it.  the method fits log(E(b)) j=Ax, (log of the expected value) not log (b)=Ax (log of the data).  

This is a critical distinction that allows direct fitting by maximum likelihood when there are zeros in the data.

Rather than an answer, I have a few questions which are relevant to finding efficiencies  for solving your problems.  You say that one of the complications is that you have many such systems.  So my first question is:  Does the coefficient matrix A vary from problem to problem?  Secondly, is the known matrix A of full column rank?  And finally, is there any additional information you have about the data  b such as error estimates or distribution?  This would be useful information in designing an algorithm.

Dear John and Thomas, thank you for the valuable comments: they have really pushed me in the right direction!

1. Thomas, regarding your questions:

The data in b is theoretically from a truncated multinomial distribution, i.e. one of the multinomial "bins" is unobservable, and so the remaining n bin counts (vector b) also have multinomial distribution but with "binomially random" total number of tries.

Matrix A is constant throughout one experiment, always with full column rank. Optimization for constant A for the solution based on QR-decomposition is obvious and was already implemented (it boils down to pre-calculating the pseudoinverse of A and multiplying it by b-s to obtain the solution). But with GLM approach (see below) the optimization is not clear for me yet.

2. John, I have extensively studied the hardcore statistical papers on capture-recapture analysis and found that the standard method for them is to formulate the models as GLMs with logarithm link function and to fit them as Poisson regression, i.e. assuming that the elements in vector b are independent and have Poisson distributions (with different lambda-parameters, of course).

As I understand it, the problem here is that the elements in b may be regarded independent "Poissons" only if all the expected occupation numbers are big enough (including that unobservable one) - and this is clearly not my case. There are publications on capture-recapture that claim that "the Poisson estimator is almost identical to the conditional multinomial estimator", but I have not read them yet (should do soon). Anyway, it would be nice to have a GLM fitting procedure that supposes multinomial distribution of the data instead of independent "Poissons". Although, because the standard GLM-fitting functions offer "Poisson" and "binomial" distributions in this class, I don't think that multinomial fitters exist...

As far as I remember from my past studies many years ago.

Your type of equations are the result of first order differential equations:

X’(t) = A X(t) with some initial condition X(0) = X0 = [X10, X20, …. Xn0]

Then we used to Diagonalize the matrix which must be an n × n matrix by the way.

Then use the equation P-1AP =λI. And do all the expansions.

Anyway after reading the theory, you will note that the exponential function of AX is only defined for square matrices AX.

Because in the expansion of the exponential eAX:

= 1 + AX + ½! (AX)2 + 1/3! (AX)3 +… = ∑∞k=0 (AX)k/k!

Is not defined for non-square matrices, because how do we define powers of AX, if A is an m × n matrix and not an n × n matrix.

Please let me know if I did help or is this a new study that I am not aware of.

Following the definition for exp^A, it exists only A is a square matrix (V. I. Arnold, ODE, M, 1971, p 92). In Your case A is a (m,n) matrix, m

Reza and Anna, I intentionally pointed out in the opening post that exp(.) here is understood as element-wise operation on a vector.

And I have already semi-solved the problem by using the GLM fitting method. It does alter the problem's formulation a bit, but it really works.

Thanks for explaining your problem more clearly. I understand your work better now, the method must be an approximation method for solutions to  that type exponential Matrix equation. I have seen the theory in some reference books, but its not my area of study. My study in this area did not go beyond the work with square matrices. After that I moved on toward computational Group theory and then topological Groups and combinatorial  group theory.