[petsc-users] Linear Least Square with regularization using petsc4py

[Romain Jolivet]

> On Jul 30, 2015, at 15: 39, Matthew Knepley <knepley at gmail.com> wrote:
> 
> On Tue, Jul 28, 2015 at 5:55 AM, Romain Jolivet <jolivetinsar at gmail.com <mailto:jolivetinsar at gmail.com>> wrote:
> Dear PETSc developers and users,
> 
> I am trying to solve a linear least square problem that requires some regularization. My plan is to use a Tikhonov scheme as a first try.
> Hence, I would like to minimise the function S(x) defined as:
> 
> S(x) = ||Ax - b||2 + alpha || Cx ||
> 
> where A is my operator (rectangular matrix with far more lines than columns), x is the vector of unknown, b is some data vector, alpha is a damping factor and C some operator of my own.
> Ideally C is a gaussian form of a covariance matrix. But at first, I would like to check what happens with an identity (good), gradient (better) or laplacian (best) operator.
> There is a mention to a MatShift function in some previous post in the mailing list concerning Tikhonov regularization, but I don’t see how to incorporate that into my python set of tools that I have set up to solve my problem (would there be a useful example somewhere?).
> 
> It would be nice to have an example for this, but I have not made it there yet. For now, it looks like
> what you want is
> 
>   A.shift(aC)
> 
> where A is the Mat that you are currently using in your linear solve, and aC is alpha*C. Does this make sense?

Thanks for the answer. 
I am not sure I get it: C is a square matrix the size of the solution vector. How does it fit in the method .shift (I thought this one needs a petsc scalar for input)?

R


>   Thanks,
> 
>      Matt
>  
> I feel like there is some possibilities using the TAO wrappers of petsc4py, but I have no idea how to get this guy to work (how do I create the TAO_objective, TAO_gradient, etc?).
> 
> To get a better understanding of my problem, my ultimate goal is to minimise the function S(x) defined as
> 
> S(x) = 1/2. * [(Ax-b)’ inv(D) (Ax-b) + (x-x0)’ inv(M) (x-x0)],
> 
> where D is some covariance matrix describing noise distribution on b, M is the covariance matrix describing the correlations between the parameters of x and x0 is the starting point.
> But that seems a bit complex so far and I will stick to the first problem yet.
> 
> Cheers,
> Romain
> 
> 
> 
> —————————————————————————————————————
> —————————————————————————————————————
> Romain Jolivet
> Postdoctoral Fellow
> 
> University of Cambridge
> Department of Earth Sciences
> Bullard Labs
> Madingley Rise
> Madingley Road
> Cambridge CB3 0EZ
> United Kingdom
> 
> email: rpj29 at cam.ac.uk <mailto:rpj29 at cam.ac.uk>
> Phone: +44 1223 748 938 <tel:%2B44%201223%20748%20938>
> Mobile: +44 7596 703 148 <tel:%2B44%207596%20703%20148>
> 
> France: +33 6 52 91 76 39 <tel:%2B33%206%2052%2091%2076%2039>
> US: +1 (626) 560 6356 <tel:%2B1%20%28626%29%20560%206356>
> —————————————————————————————————————
> —————————————————————————————————————
> 
> 
> 
> 
> -- 
> What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.
> -- Norbert Wiener

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20150731/06709b5b/attachment.html>