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

Romain Jolivet jolivetinsar at gmail.com
Tue Jul 28 05:55:40 CDT 2015

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


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

More information about the petsc-users mailing list