What you have here is a "blind deconvolution" problem. There are strategies for this, but it's a very difficult problem and depends on having some useful prior information.

Are these draws independent, or correlated in time? If they're correlated, do you know anything about the correlation?

Let's be clear about "all I know is the distribution of $\eta$." Are you saying that you know a probability distribution for $\eta(t_{i})$ at each point in time $t_{i}$, or that you also know the correlation between $\eta(t_{i})$ at different points in time?

What is the probabilistic interpretation of the entries of $A^{-1}$? Are you looking for the entries of a covariance matrix of the form $(X^{T}X)^{-1}$ for example? "Almost Block Diagonal" really isn't gonig to be enough to help here- as soon as you've got one entry linking two blocks, you can get fill in of the linking block in the inverse matrix.

The point is that some reviewer might well ask to have the paper rewritten using the common convention and it's quite likely that an editor would agree that the author should do this. If this happens it will slow down the process of getting your paper accepted.

The operator $O$ here is computed by a numerical simulation that may or may not be noisy. If the simulation results are noisy, then using finite difference derivatives may be a poor choice. An alternative to consider in that case is the use of response surface methods. However, if the simulation results are reasonably good then finite difference derivatives are likely to work just fine.

Gradient descent converges at only a linear rate, whereas Newton's method (and LM is essentially Newton's method) converges at a quadratic rate. LM can be much faster in practice.

One important issue with the LM method is that it a local search method. In some circumstances, if your least squares problem has local minima it may return a solution that is locally but not globally optimal. It's possible in that case that you'd miss a solution that really did satisfy $O(x,y)\psi=\phi$. A simple approach that can help with this problem is to start LM running from many different initial guesses and then take the best solution found.

The set of full rank matrices is not convex, and I'm not aware of any work on optimizing over the set of full rank matrices. I think you're better off staying with the original formulation. For reasonably small values of $n$, the original formulation can readily be attacked by branch and bound methods. How big is your $n$?

Let me also mention that if $\alpha$ is fixed and you want to vary $s$, then you might find that a better way to go is to compute the eigenvalue decomposition of $Q-\alpha J$, and then adjust for $sI$ by adding $s$ to the eigenvalues.

A classic reference (and the paper is available online as a free .pdf) is: P. Gill, G. Golub, W. Murray, and M. Saunders. Methods for modifying matrix factorizations. Mathematics of Computation, 126(28):505-535, 1974. stanford.edu/group/SOL/papers/ggms74.pdf