Return-Path: X-Original-To: apmail-commons-dev-archive@www.apache.org Delivered-To: apmail-commons-dev-archive@www.apache.org Received: from mail.apache.org (hermes.apache.org [140.211.11.3]) by minotaur.apache.org (Postfix) with SMTP id EA47A74D3 for ; Sun, 7 Aug 2011 16:46:04 +0000 (UTC) Received: (qmail 23043 invoked by uid 500); 7 Aug 2011 16:46:04 -0000 Delivered-To: apmail-commons-dev-archive@commons.apache.org Received: (qmail 22701 invoked by uid 500); 7 Aug 2011 16:46:03 -0000 Mailing-List: contact dev-help@commons.apache.org; run by ezmlm Precedence: bulk List-Help: List-Unsubscribe: List-Post: List-Id: Reply-To: "Commons Developers List" Delivered-To: mailing list dev@commons.apache.org Received: (qmail 22693 invoked by uid 99); 7 Aug 2011 16:46:03 -0000 Received: from nike.apache.org (HELO nike.apache.org) (192.87.106.230) by apache.org (qpsmtpd/0.29) with ESMTP; Sun, 07 Aug 2011 16:46:03 +0000 X-ASF-Spam-Status: No, hits=1.5 required=5.0 tests=FREEMAIL_FROM,HTML_MESSAGE,RCVD_IN_DNSWL_LOW,SPF_PASS X-Spam-Check-By: apache.org Received-SPF: pass (nike.apache.org: domain of gsterijevski@gmail.com designates 209.85.220.171 as permitted sender) Received: from [209.85.220.171] (HELO mail-vx0-f171.google.com) (209.85.220.171) by apache.org (qpsmtpd/0.29) with ESMTP; Sun, 07 Aug 2011 16:45:54 +0000 Received: by vxh13 with SMTP id 13so1921288vxh.30 for ; Sun, 07 Aug 2011 09:45:33 -0700 (PDT) DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=mime-version:in-reply-to:references:date:message-id:subject:from:to :content-type; bh=pd8PFszcjY4CQCeNDCmn5UGWwPyRA+v7QkOOFlUo2xk=; b=K7poM0tx7Re18RMdClcN114HvAkikwh7Qz5sXcH4ptiqSjSClhMrrhMn7Il7psVwnT ftGkjYc/coFTHXlwFND+LYxOTBvL5mTadMKE7wCZ++j+gciIWl3kyMTpR/S5UQkGHoJA k5QQ0joxXlWHCefqulCyfO4u5QuuWea5kzBXM= MIME-Version: 1.0 Received: by 10.220.8.193 with SMTP id i1mr1333225vci.108.1312735533647; Sun, 07 Aug 2011 09:45:33 -0700 (PDT) Received: by 10.220.95.136 with HTTP; Sun, 7 Aug 2011 09:45:33 -0700 (PDT) In-Reply-To: References: <4E27BBD9.9060805@gmail.com> <4E2846AF.3060904@gmail.com> <4E298977.4030708@gmail.com> <4E2AF6F6.8060301@gmail.com> <4E2B37EB.6050204@gmail.com> Date: Sun, 7 Aug 2011 11:45:33 -0500 Message-ID: Subject: Re: [math] Pivoting in QR decomposition From: Greg Sterijevski To: Commons Developers List Content-Type: multipart/alternative; boundary=bcaec54ee99e2ac6dd04a9ed0f6b X-Virus-Checked: Checked by ClamAV on apache.org --bcaec54ee99e2ac6dd04a9ed0f6b Content-Type: text/plain; charset=ISO-8859-1 Your selection methodology might be more complicated than pairwise comparison. So I think you would have to pass the entire set of remaining columns. That is an implementation detail, however. On the interface, yes the public face should not expose the constructor which takes the PivotStrategyComparator. It should be part of a protected constructor or, I suppose, we could have a static factory class which spits out different implementations. -Greg On Sun, Aug 7, 2011 at 11:33 AM, Chris Nix wrote: > To avoid having the strategy method being overriden in the constructor, we > could instead implement a constructor that takes a column Comparator that > determines if two columns should be exchanged at each stage. > > In the interest of maintaining a clean interface, I don't know if this > constructor should be public though. > > Chris > > On 7 August 2011 16:43, Greg Sterijevski wrote: > > > Chris, > > > > In regard to the pivoting, do you think that it would be useful to allow > > subclasses to pivot using other strategies? The pivoting I have looks for > > the next column with the largest norm. For most garden variety problems > > this > > will be okay. However, you can (I am not sure how effective this will be) > > choose the column which has the lowest [average] inner product (ie is > least > > correlated in some sense to the remaining columns). > > > > The easiest way to accomplish this is to present the remainder columns to > > some method, and have that method return an index... The method would be > > protected so that it could be overridden by classes wishing to modify > that > > behavior. The only problem I see is that the actual decomposition would > > need > > to be moved out of the constructor (we would have an overridable method > > being called in the constructor). > > > > -Greg > > > > On Sun, Aug 7, 2011 at 7:29 AM, Chris Nix wrote: > > > > > Oooops, the ********* below should read MATH-630. > > > > > > Chris > > > > > > On 7 August 2011 13:28, Chris Nix wrote: > > > > > > > Thanks, Greg, for looking more at this. Apologies I've not been able > > to > > > > focus on this too much recently. > > > > > > > > Yes, the approach above works, however the solvers also require a > > change > > > so > > > > that they don't unexpectedly solve for a permuted matrix. > > Additionally, > > > a > > > > getRank() method can now be implemented much more efficiently than > the > > > > getRank() from SingularValueDecomposition. > > > > > > > > I submitted an initial patch at *********** with these bits working, > > > > however this patch introduces a new RRQRDecomposition interface > > extending > > > > QRDecomposition. In light of the insights above from the team, I'll > > > > implement it instead within the existing class structure and > re-submit. > > > > > > > > If the pivoting code is to be included in QRDecomposition, then > perhaps > > > an > > > > extra constructor is required to perform column pivoting since the > RRQR > > > > decomposition of matrix A produces Q, R and P such A = QRP^T and > would > > > break > > > > any existing code that expects a plain decomposition of A=QR. > > > > > > > > Chris. > > > > > > > > On 6 August 2011 21:34, Greg Sterijevski > > wrote: > > > > > > > >> Hello All, > > > >> > > > >> Not sure if this is stepping on toes, but here is what I thought of > to > > > >> deal > > > >> with pivoting. I would only need to modify the constructor of the > > > >> QRDecompImpl class: > > > >> > > > >> public QRDecompositionPivotImpl(RealMatrix matrix) { > > > >> > > > >> final int m = matrix.getRowDimension(); > > > >> final int n = matrix.getColumnDimension(); > > > >> final int mn = FastMath.min(m, n); > > > >> > > > >> qrt = matrix.transpose().getData(); > > > >> rDiag = new double[FastMath.min(m, n)]; > > > >> cachedQ = null; > > > >> cachedQT = null; > > > >> cachedR = null; > > > >> cachedH = null; > > > >> > > > >> > > > >> /* > > > >> * The QR decomposition of a matrix A is calculated using > > > >> Householder > > > >> * reflectors by repeating the following operations to each > > minor > > > >> * A(minor,minor) of A: > > > >> */ > > > >> > > > >> pivots = new int[n]; > > > >> for (int i = 0; i < n; i++) { > > > >> pivots[i] = i; > > > >> } > > > >> > > > >> int pivotIdx = -1; > > > >> double bestNorm = 0.0; > > > >> for (int minor = 0; minor < mn; minor++) { > > > >> bestNorm = 0.0; > > > >> pivotIdx = 0; > > > >> for (int i = minor; i < n; i++) { > > > >> final double[] qrtMinor = qrt[i]; > > > >> double xNormSqr = 0.0; > > > >> for (int row = minor; row < m; row++) { > > > >> final double c = qrtMinor[row]; > > > >> xNormSqr += c * c; > > > >> } > > > >> if (xNormSqr > bestNorm) { > > > >> bestNorm = xNormSqr; > > > >> pivotIdx = i; > > > >> } > > > >> } > > > >> > > > >> > > > >> if( pivotIdx != minor){ > > > >> int pvt = pivots[minor]; > > > >> pivots[minor] = pivots[pivotIdx]; > > > >> double[] qswp = qrt[minor]; > > > >> qrt[minor] = qrt[pivotIdx]; > > > >> for( int i = minor + 1; i < n; i++){ > > > >> if( i <= pivotIdx ){ > > > >> final int itmp = pivots[i]; > > > >> pivots[i] = pvt; > > > >> pvt = itmp; > > > >> final double[] tmp = qrt[i]; > > > >> qrt[i] = qswp; > > > >> qswp=tmp; > > > >> } > > > >> } > > > >> } > > > >> > > > >> final double[] qrtMinor = qrt[minor]; > > > >> /* > > > >> * Let x be the first column of the minor, and a^2 = > |x|^2. > > > >> * x will be in the positions qr[minor][minor] through > > > >> qr[m][minor]. > > > >> * The first column of the transformed minor will be > > > >> (a,0,0,..)' > > > >> * The sign of a is chosen to be opposite to the sign of > > the > > > >> first > > > >> * component of x. Let's find a: > > > >> */ > > > >> final double a = (qrtMinor[minor] > 0) ? > > > >> -FastMath.sqrt(bestNorm) : FastMath.sqrt(bestNorm); > > > >> rDiag[minor] = a; > > > >> > > > >> if (a != 0.0) { > > > >> > > > >> /* > > > >> * Calculate the normalized reflection vector v and > > > >> transform > > > >> * the first column. We know the norm of v > beforehand: > > v > > > = > > > >> x-ae > > > >> * so |v|^2 = = -2a+a^2 = > > > >> * a^2+a^2-2a = 2a*(a - ). > > > >> * Here is now qr[minor][minor]. > > > >> * v = x-ae is stored in the column at qr: > > > >> */ > > > >> qrtMinor[minor] -= a; // now |v|^2 = > > > -2a*(qr[minor][minor]) > > > >> > > > >> /* > > > >> * Transform the rest of the columns of the minor: > > > >> * They will be transformed by the matrix H = > > > I-2vv'/|v|^2. > > > >> * If x is a column vector of the minor, then > > > >> * Hx = (I-2vv'/|v|^2)x = x-2vv'x/|v|^2 = x - > > > 2/|v|^2 > > > >> v. > > > >> * Therefore the transformation is easily calculated > by > > > >> * subtracting the column vector (2/|v|^2)v from > > x. > > > >> * > > > >> * Let 2/|v|^2 = alpha. From above we have > > > >> * |v|^2 = -2a*(qr[minor][minor]), so > > > >> * alpha = -/(a*qr[minor][minor]) > > > >> */ > > > >> for (int col = minor + 1; col < n; col++) { > > > >> final double[] qrtCol = qrt[col]; > > > >> double alpha = 0; > > > >> for (int row = minor; row < m; row++) { > > > >> alpha -= qrtCol[row] * qrtMinor[row]; > > > >> } > > > >> alpha /= a * qrtMinor[minor]; > > > >> > > > >> // Subtract the column vector alpha*v from x. > > > >> for (int row = minor; row < m; row++) { > > > >> qrtCol[row] -= alpha * qrtMinor[row]; > > > >> } > > > >> } > > > >> } > > > >> } > > > >> } > > > >> > > > >> > > > >> All I am doing is looking forward and taking the next column with > the > > > >> largest norm. Then I rearrange the Q's. Is this what you had in mind > > > >> Chris? > > > >> The result will be Q,R and the pivot array for which we can > implement > > a > > > >> getter? > > > >> > > > >> -Greg > > > >> > > > > > > > > > > > > > > --bcaec54ee99e2ac6dd04a9ed0f6b--