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

[Matthew Knepley]

On Fri, Jul 31, 2015 at 3:30 AM, Romain Jolivet <jolivetinsar at gmail.com>
wrote:

>
> 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>
> 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)?
>

Okay, could you then write down the linear algebra you want to do? Then I
can tell you which operations to use.

  Thanks,

     Matt


> 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
>> Phone: +44 1223 748 938
>> Mobile: +44 7596 703 148
>>
>> France: +33 6 52 91 76 39
>> US: +1 (626) 560 6356
>> —————————————————————————————————————
>> —————————————————————————————————————
>>
>>
>
>
> --
> 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
>
>
>


-- 
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/d14fbccd/attachment.html>