Re: [petsc-users] L1 or projection type regularization with PETSc
Lingyun Qiu writes: > On Thu, Apr 7, 2016 at 2:24 PM, Jed Brown wrote: > >> Lingyun Qiu writes: >> >> > Dear all, >> > >> > I am working an optimization problem as >> > min_x ||Ax - b||_2^2 + alpha ||x||_1 >> > >> > For the fidelity term, we use L2 norm. We use L1 norm for the >> > regularization term. Without regularization term, i.e., alpha=0, we >> > iteratively solve the problem as >> > x_k+1 = KSP(x_k). >> >> I.e., a single KSPSolve because the optimization problem is quadratic? > > Yes. A single KSPSolve because of the quadratic form. k stands for the > iteration index inside KSPSolve. But the k+1 iterate is not a function of just the previous iterate, it is a function of all previous iterates (sometimes via recurrence). >> > I plan to use the split Bregman method to solve the regularized problem. >> It >> > reads as, >> > y_k+1 = KSP(x_k) >> > x_k+1 = B(y_k+1) >> > Here B() is the function related to the Bregman method. It works as a >> > post-processing of the iterates. >> >> But outside the KSPSolve, not at each iteration of the KSPSolve... >> > I need to modify the iterate after each iteration inside the KSPSolve. Yes, > B is a nonlinear operation. That nonlinear operation breaks orthogonality, so your Krylov method won't work. You can't just tinker willy-nilly with the vectors inside the Krylov solver. There are specialized Krylov methods for solving problems with certain constraints, such as QCG for solving SPD systems subject to trust region constraints. But not in general. Call KSPSolve and then do your shrink operation, rinse and repeat, perhaps with nonlinear acceleration. signature.asc Description: PGP signature
Re: [petsc-users] L1 or projection type regularization with PETSc
I think *KSPSetPostSolve * is to add a post-processing after all the
iterations inside KSPSolve finish. I need to modify the iterates inside
KSPSolve or the linear solver.
I just realize that this nonlinear operation will also affect the Krylov
subspaces. Perhaps I should look for some other linear solver instead of
KSPSolve. Comments are welcome.
On Thu, Apr 7, 2016 at 2:24 PM, Matthew Knepley wrote:
> On Thu, Apr 7, 2016 at 2:13 PM, Lingyun Qiu wrote:
>
>> Dear all,
>>
>> I am working an optimization problem as
>> min_x ||Ax - b||_2^2 + alpha ||x||_1
>>
>> For the fidelity term, we use L2 norm. We use L1 norm for the
>> regularization term. Without regularization term, i.e., alpha=0, we
>> iteratively solve the problem as
>> x_k+1 = KSP(x_k).
>>
>> I plan to use the split Bregman method to solve the regularized problem.
>> It reads as,
>> y_k+1 = KSP(x_k)
>> x_k+1 = B(y_k+1)
>> Here B() is the function related to the Bregman method. It works as a
>> post-processing of the iterates.
>>
>> I am wondering is there a way to combine this post-processing with the
>> KSP solver? A brute-force way is modify the initial guess and set the max
>> iteration number to 1.
>>
>
> Are you asking for something like this:
>
>
> http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/KSP/KSPSetPostSolve.html
>
> Thanks,
>
> Matt
>
>
>> This is also related to the projection type regularization:
>> min_{x in subspace G} ||Ax-b||^2_2
>> The scheme is
>> y_k+1 = KSP(x_k)
>> x_k+1 = P_G(y_k+1)
>> where P_G is the projection to subspace G.
>>
>> Lingyun Qiu
>>
>
>
>
> --
> 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
>
Re: [petsc-users] L1 or projection type regularization with PETSc
On Thu, Apr 7, 2016 at 2:24 PM, Jed Brown wrote: > Lingyun Qiu writes: > > > Dear all, > > > > I am working an optimization problem as > > min_x ||Ax - b||_2^2 + alpha ||x||_1 > > > > For the fidelity term, we use L2 norm. We use L1 norm for the > > regularization term. Without regularization term, i.e., alpha=0, we > > iteratively solve the problem as > > x_k+1 = KSP(x_k). > > I.e., a single KSPSolve because the optimization problem is quadratic? Yes. A single KSPSolve because of the quadratic form. k stands for the iteration index inside KSPSolve. > > > I plan to use the split Bregman method to solve the regularized problem. > It > > reads as, > > y_k+1 = KSP(x_k) > > x_k+1 = B(y_k+1) > > Here B() is the function related to the Bregman method. It works as a > > post-processing of the iterates. > > But outside the KSPSolve, not at each iteration of the KSPSolve... > I need to modify the iterate after each iteration inside the KSPSolve. Yes, B is a nonlinear operation. > > B is a nonlinear operation, so you can't put it inside KSP. You could > put it inside SNES (perhaps with nonlinear preconditioning, NGMRES or > QN), but I don't think that's what you're asking for here. > I will check these functions. Thanks for the quick reply.
Re: [petsc-users] L1 or projection type regularization with PETSc
On Thu, Apr 7, 2016 at 2:13 PM, Lingyun Qiu wrote:
> Dear all,
>
> I am working an optimization problem as
> min_x ||Ax - b||_2^2 + alpha ||x||_1
>
> For the fidelity term, we use L2 norm. We use L1 norm for the
> regularization term. Without regularization term, i.e., alpha=0, we
> iteratively solve the problem as
> x_k+1 = KSP(x_k).
>
> I plan to use the split Bregman method to solve the regularized problem.
> It reads as,
> y_k+1 = KSP(x_k)
> x_k+1 = B(y_k+1)
> Here B() is the function related to the Bregman method. It works as a
> post-processing of the iterates.
>
> I am wondering is there a way to combine this post-processing with the KSP
> solver? A brute-force way is modify the initial guess and set the max
> iteration number to 1.
>
Are you asking for something like this:
http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/KSP/KSPSetPostSolve.html
Thanks,
Matt
> This is also related to the projection type regularization:
> min_{x in subspace G} ||Ax-b||^2_2
> The scheme is
> y_k+1 = KSP(x_k)
> x_k+1 = P_G(y_k+1)
> where P_G is the projection to subspace G.
>
> Lingyun Qiu
>
--
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
Re: [petsc-users] L1 or projection type regularization with PETSc
Lingyun Qiu writes: > Dear all, > > I am working an optimization problem as > min_x ||Ax - b||_2^2 + alpha ||x||_1 > > For the fidelity term, we use L2 norm. We use L1 norm for the > regularization term. Without regularization term, i.e., alpha=0, we > iteratively solve the problem as > x_k+1 = KSP(x_k). I.e., a single KSPSolve because the optimization problem is quadratic? > I plan to use the split Bregman method to solve the regularized problem. It > reads as, > y_k+1 = KSP(x_k) > x_k+1 = B(y_k+1) > Here B() is the function related to the Bregman method. It works as a > post-processing of the iterates. But outside the KSPSolve, not at each iteration of the KSPSolve... B is a nonlinear operation, so you can't put it inside KSP. You could put it inside SNES (perhaps with nonlinear preconditioning, NGMRES or QN), but I don't think that's what you're asking for here. signature.asc Description: PGP signature
[petsc-users] L1 or projection type regularization with PETSc
Dear all,
I am working an optimization problem as
min_x ||Ax - b||_2^2 + alpha ||x||_1
For the fidelity term, we use L2 norm. We use L1 norm for the
regularization term. Without regularization term, i.e., alpha=0, we
iteratively solve the problem as
x_k+1 = KSP(x_k).
I plan to use the split Bregman method to solve the regularized problem. It
reads as,
y_k+1 = KSP(x_k)
x_k+1 = B(y_k+1)
Here B() is the function related to the Bregman method. It works as a
post-processing of the iterates.
I am wondering is there a way to combine this post-processing with the KSP
solver? A brute-force way is modify the initial guess and set the max
iteration number to 1.
This is also related to the projection type regularization:
min_{x in subspace G} ||Ax-b||^2_2
The scheme is
y_k+1 = KSP(x_k)
x_k+1 = P_G(y_k+1)
where P_G is the projection to subspace G.
Lingyun Qiu
[petsc-users] PDESoft 2016, UK, 4-8 Jul 2016 (abstract submission deadline extended)
The abstract submission deadline has been extended to April 22. Please consider this small and interesting conference on your summer schedule. We are happy to announce that PDESoft 2016 will be held at the University of Warwick, UK. This main conference will run from 4-6 July and the remaining time will be dedicated to small group coding days. PDESoft 2016 provides a forum for the developers of open-source tools solving the diverse stages of the numerical PDE process to exchange the latest developments and future ideas. If you develop: - meshing tools; - HPC tools; - solvers for large systems of equations; - numerical PDE solvers; - visualisation systems for PDE; - coupling with or user interfaces to scientific software; or any other part of the PDE toolchain, this is the workshop for you. For more details please visit http://warwick.ac.uk/pdesoft2016 signature.asc Description: PGP signature
Re: [petsc-users] Ordering with MatCreateNest and MPI
Alessandro Francavilla writes: > Thanks for the reply; probably I've not been clear enough. I don't need to > separate out components, I just need to understand the underlying orderings. > > For instance, suppose that I want to multiply the matrix with a null vector > with a 1 in the i-th location (the result is basically the i-th column of > the matrix): > > CALL > MatCreateNest(PETSC_COMM_WORLD,2,PETSC_NULL_OBJECT,2,PETSC_NULL_OBJECT,matAr > ray,A,IERR) > CALL MatCreateVecs(A,x,b,IERR) > CALL VecSet(x,0.+0.*PETSC_i,IERR) > CALL VecSetValues(x,1,i-1,(1.+0.*PETSC_i),INSERT_VALUES,IERR) > CALL VecAssemblyBegin(x,IERR) > CALL VecAssemblyEnd(x,IERR) > CALL MatMult(A,x,b,IERR) > > The above code works perfectly if I execute it single process (mpirun -np > 1); if I run it on >1 MPI processes, the result is not what I would expect: > rather than returning the i-th column, it returns (e.g., if I call VecView > on b) some permutation of a different column (say a permutation of column > j). > > Also, if I substitute the call to MatCreateNest() with > A = matArray(1) But this is a different matrix. Why would you expect that the MatNest has the same indices as its first Mat? > Then everything works as expected, for arbitrary numbers of processes. By > working as expected I mean that calling VecView on b prints the i-th column > of A with its natural ordering. > > So, my question is: > 1) how do I generate the input vector x such that the MatMult operation > returns the i-th column of A? It is returning the ith column of A. A is not just matArray(1). > 2) how do I reorder the result of MatMult such that I get the i-th column > with its "natural" ordering? > > Maybe I'm even misunderstanding the output of VecView, which prints the > elements process by process, and thus the order corresponds to the natural > ordering only if the elements are distributed consecutively among processes? > Which would explain why I'm seeing a permutation of the column. However, I'm > seeing the permutation of a DIFFERENT column of the matrix, which means that > the input vector does not have a 1 in the i-th position when I fill it with > VecSetValues(x,1,i-1,(1.+0.*PETSC_i),INSERT_VALUES,IERR). PETSc Vecs are distributed with contiguous blocks assigned to each process. Mat's column and row indices are induced by their action on Vecs. MatNest does not move matrix data around, so its blocks produce the portions of Vecs that are owned on the current process. Consequently, MatNest induces an ordering in which one field gets the first chunk of a contiguous space on each process, with the intervening entries for the other field(s). If you have two matrices, B and C, each distributed over the same two processes, then the MatNest [B 0; 0 C] will keep both B and C distributed. I.e., it does not move all of B to rank 0 and all of C to rank 1, or some funny business like that. If you don't want to think about this ordering, extract the subvector the way I described in the previous email. signature.asc Description: PGP signature
Re: [petsc-users] Ordering with MatCreateNest and MPI
Thanks for the reply; probably I've not been clear enough. I don't need to separate out components, I just need to understand the underlying orderings. For instance, suppose that I want to multiply the matrix with a null vector with a 1 in the i-th location (the result is basically the i-th column of the matrix): CALL MatCreateNest(PETSC_COMM_WORLD,2,PETSC_NULL_OBJECT,2,PETSC_NULL_OBJECT,matAr ray,A,IERR) CALL MatCreateVecs(A,x,b,IERR) CALL VecSet(x,0.+0.*PETSC_i,IERR) CALL VecSetValues(x,1,i-1,(1.+0.*PETSC_i),INSERT_VALUES,IERR) CALL VecAssemblyBegin(x,IERR) CALL VecAssemblyEnd(x,IERR) CALL MatMult(A,x,b,IERR) The above code works perfectly if I execute it single process (mpirun -np 1); if I run it on >1 MPI processes, the result is not what I would expect: rather than returning the i-th column, it returns (e.g., if I call VecView on b) some permutation of a different column (say a permutation of column j). Also, if I substitute the call to MatCreateNest() with A = matArray(1) Then everything works as expected, for arbitrary numbers of processes. By working as expected I mean that calling VecView on b prints the i-th column of A with its natural ordering. So, my question is: 1) how do I generate the input vector x such that the MatMult operation returns the i-th column of A? 2) how do I reorder the result of MatMult such that I get the i-th column with its "natural" ordering? Maybe I'm even misunderstanding the output of VecView, which prints the elements process by process, and thus the order corresponds to the natural ordering only if the elements are distributed consecutively among processes? Which would explain why I'm seeing a permutation of the column. However, I'm seeing the permutation of a DIFFERENT column of the matrix, which means that the input vector does not have a 1 in the i-th position when I fill it with VecSetValues(x,1,i-1,(1.+0.*PETSC_i),INSERT_VALUES,IERR). Thanks, Alessandro -Original Message- From: Jed Brown [mailto:j...@jedbrown.org] Sent: giovedì 7 aprile 2016 15:26 To: FRANCAVILLA MATTEO ALESSANDRO; petsc-users@mcs.anl.gov Subject: Re: [petsc-users] Ordering with MatCreateNest and MPI FRANCAVILLA MATTEO ALESSANDRO writes: > Hi, > > I'm trying to use a matrix composed of different sub-matrix blocks; I > create the matrix with MatCreateNest(). I'm finding it hard to > understand the ordering (in MPI applications) of the result of MatMult operations. > Apparently the result of the multiplication is correct, except that > both input and output vectors go through some sort of permutation. How > should I use the MatMult routine, i.e. how to pass the input vector > with the correct ordering and how to get the correct ordering of the result? > > Below a simple Fortran example to demonstrate what I'm trying to do > (running the example with 1 process works fine, with 2 or more > processes I don't understand how to get the correct ordering). When > using standard matrices (i.e., not MatNest) I've never had to care about orderings. MatNest is semantically identical to other matrix formats. In fact, we strongly recommend that you write your code so that your assembly and all subsequent operations do not depend in any way on the matrix being MATNEST of some other format. See src/snes/examples/tutorials/ex28.c for an example of this. But anyway, you can use VecGetSubVector to separate out the components (MatNestGetISs to get the index sets if you don't already have them).
Re: [petsc-users] Prelloc: Get coordinates of local diagonal submatrix
Florian Lindner writes: >Am I right to assume that nnz is to in the the PETSc ordering of rows and not >application ordering? Yes signature.asc Description: PGP signature
Re: [petsc-users] Prelloc: Get coordinates of local diagonal submatrix
Hello,
> MatGetOwnershipRangesColumn
Thanks! That's what I was looking for.
I have set a MatSetLocalToGlobalMapping on the matrix and I use
MatMPIAIJSetPreallocation(_matrixA.matrix, 0, d_nnz, 0, o_nnz)
to set the preallocation. Am I right to assume that nnz is to in the the PETSc
ordering of rows and not application ordering?
Best Thanks,
Florian
Am Mittwoch, 6. April 2016, 10:19:23 CEST schrieben Sie:
> On Wed, Apr 6, 2016 at 10:12 AM, Florian Lindner
> wrote:
>
> > Hello,
> >
> > in order to preallocate I have to know whether a non-zero will be in the
> > local diagonal submatrix or not.
> >
> >
> > http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatMPIAIJSetPreallocation.html
> >
> > says:
> >
> > The DIAGONAL portion of the local submatrix of a processor can be defined
> > as the submatrix which is obtained by extraction the part corresponding to
> > the rows r1-r2 and columns c1-c2 of the global matrix, where r1 is the
> > first row that belongs to the processor, r2 is the last row belonging to
> > the this processor, and c1-c2 is range of indices of the local part of a
> > vector suitable for applying the matrix to. This is an mxn matrix. In the
> > common case of a square matrix, the row and column ranges are the same and
> > the DIAGONAL part is also square. The remaining portion of the local
> > submatrix (mxN) constitute the OFF-DIAGONAL portion.
> >
> > Therefore I compute the begin and end rows and cols like that:
> >
> > int ownerRangeBeginA = _matrixA.ownerRange().first;
> > int ownerRangeEndA = _matrixA.ownerRange().second;
> >
> > // Uses: MatGetVecs(_matrixA.matrix, &vector, nullptr);
> > petsc::Vector diagTest{_matrixA, "test", petsc::Vector::RIGHT};
> >
> > int localDiagColBegin = diagTest.ownerRange().first;
> > int localDiagColEnd = diagTest.ownerRange().second;
> >
> > Debug("Local Submatrix Rows = " << ownerRangeBeginA << " / " <<
> > ownerRangeEndA <<
> > ", Local Submatrix Cols = " << localDiagColBegin << " / " <<
> > localDiagColEnd);
> >
> >
> > It's a little bit tainted by my PETSc C++ helper function, but I think you
> > get the code.
> >
> > Is there a way to do it more elegantly? Without instantiating a vector
> > just for the purpose of gettings its owner range? And then testing a if an
> > index is inside the local diagonal submatrix?
> >
>
> http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatGetOwnershipRangesColumn.html
>
> Thanks,
>
> Matt
>
>
> > Thanks,
> > Florian
> >
>
>
>
>
Re: [petsc-users] Ordering with MatCreateNest and MPI
FRANCAVILLA MATTEO ALESSANDRO writes: > Hi, > > I'm trying to use a matrix composed of different sub-matrix blocks; I create > the matrix with MatCreateNest(). I'm finding it hard to understand the > ordering (in MPI applications) of the result of MatMult operations. > Apparently the result of the multiplication is correct, except that both > input and output vectors go through some sort of permutation. How should I > use the MatMult routine, i.e. how to pass the input vector with the correct > ordering and how to get the correct ordering of the result? > > Below a simple Fortran example to demonstrate what I'm trying to do (running > the example with 1 process works fine, with 2 or more processes I don't > understand how to get the correct ordering). When using standard matrices > (i.e., not MatNest) I've never had to care about orderings. MatNest is semantically identical to other matrix formats. In fact, we strongly recommend that you write your code so that your assembly and all subsequent operations do not depend in any way on the matrix being MATNEST of some other format. See src/snes/examples/tutorials/ex28.c for an example of this. But anyway, you can use VecGetSubVector to separate out the components (MatNestGetISs to get the index sets if you don't already have them). signature.asc Description: PGP signature
[petsc-users] Ordering with MatCreateNest and MPI
Hi, I'm trying to use a matrix composed of different sub-matrix blocks; I create the matrix with MatCreateNest(). I'm finding it hard to understand the ordering (in MPI applications) of the result of MatMult operations. Apparently the result of the multiplication is correct, except that both input and output vectors go through some sort of permutation. How should I use the MatMult routine, i.e. how to pass the input vector with the correct ordering and how to get the correct ordering of the result? Below a simple Fortran example to demonstrate what I'm trying to do (running the example with 1 process works fine, with 2 or more processes I don't understand how to get the correct ordering). When using standard matrices (i.e., not MatNest) I've never had to care about orderings. Thanks, Alessandro PROGRAM Main !. implicit none #include #include #include !* ! PETSc variables !* Vec :: x, b Mat :: matArray(4), A PetscInt:: II, JJ, FromRow, ToRow, IMAT PetscErrorCode :: IERR PetscInt, ALLOCATABLE :: INDICES(:) PetscScalar, ALLOCATABLE :: VALUES(:) PetscInt, PARAMETER :: N = 2! size of sub-matrices CALL PetscInitialize(PETSC_NULL_CHARACTER,IERR) ALLOCATE(INDICES(N),VALUES(N)) INDICES = (/ (ii, ii=0, N-1) /) ! Form 4 blocks of size NxN DO IMAT = 1, 4 CALL MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,N,N,PETSC_NULL_SCALAR,matArray(IMAT),IERR) CALL MatGetOwnershipRange(matArray(IMAT),FromRow,ToRow,IERR) DO II = FromRow, ToRow-1 DO JJ = 1, N VALUES(JJ) = (IMAT-1)*100 + N*II + JJ END DO CALL MatSetValues(matArray(IMAT),1,II,N,INDICES,VALUES,INSERT_VALUES,IERR) END DO CALL MatAssemblyBegin(matArray(IMAT),MAT_FINAL_ASSEMBLY,IERR) CALL MatAssemblyEnd(matArray(IMAT),MAT_FINAL_ASSEMBLY,IERR) END DO ! Form a 2x2 block matrix made as [ matArray(1) matArray(2) ; matArray(3) matArray(4) ] ! If N = 2 the matrix is: ! A =[ 1 2 | 101 102 ; ! 3 4 | 102 103 ; ! ___ ! 201 202 | 401 402 ; !203 204 | 403 404 ] CALL MatCreateNest(PETSC_COMM_WORLD,2,PETSC_NULL_OBJECT,2,PETSC_NULL_OBJECT,matArray,A,IERR) CALL MatCreateVecs(A,x,b,IERR) CALL VecSet(x,0.+0.*PETSC_i,IERR) CALL VecSetValues(x,1,2,(1.+0.*PETSC_i),INSERT_VALUES,IERR) CALL VecAssemblyBegin(x,IERR) CALL VecAssemblyEnd(x,IERR) CALL VecView(x,PETSC_VIEWER_STDOUT_WORLD,IERR) CALL MatMult(A,x,b,IERR) ! A*x extracts the 3rd column of A CALL VecView(b,PETSC_VIEWER_STDOUT_WORLD,IERR) DEALLOCATE(INDICES,VALUES) CALL PetscFinalize(IERR) END PROGRAM
