linearproj, affineproj | R Documentation |
Computes the projection of points in the columns of B onto the linear subspace spaned by the columns of A, resp. the projection of a point onto an affine subspace and its distance.
linearproj(A, B)
affineproj(x0, C, b, unbound = TRUE, maxniter = 100)
A |
Matrix whose columns span a subspace of some R^n. |
B |
Matrix whose columns are to be projected. |
x0 |
Point in R^n to be projected onto C x = b. |
C , b |
Matrix and vector, defining an affine subspace as C x = b |
unbound |
Logical; require all x >= 0 if unbound is false. |
maxniter |
Maximum number of iterations (if is unbound is false). |
linearproj
projects points onto a linear subspace in R^n.
The columns of A are assumed be the basis of a linear subspace, esp.
they are required to be linearly independent. The columns of matrix B
define points in R^n that will be projected onto A, and their resp.
coefficients in terms of the basis in A are computed.
The columns of A need to be linearly independent; if not, generate an
orthonormal basis of this subspace with orth(A)
. If you want to
project points onto a subspace that is defined by A x = 0
, then
generate an orthonormal basis of the nullspace of A with null(A)
.
Technically, the orthogonal projection can be determined by a finite 'Fourier expansion' with coefficients calculated as scalar products, see the examples.
affineproj
projects (single) points onto an affine subspace
defined by A x = b
and calculates the distance of x0
from
this subspace. The calculation is based on the following formula:
p = (I - A' (A A')^{-1}) x0 + A' (A A')^{-1} b
Technically, if a
is one solution of C x = b
, then the
projection onto C can be derived from the projection onto
S = {C x = 0}
with proj_C(x) = a + proj_S(x - a)
,
see the examples.
In case the user requests the coordinates of the projected point to be positive, an iteration procedure is started where negative coordinates are set to zero in each iteration.
The functions linearproj
returns a list with components P and Q.
The columns of P contain the coefficients – in the basis of A – of the
corresponding projected points in B, and the columns of Q are the the
coordinates of these points in the natural coordinate system of R^n.
affineproj
returns a list with components proj
, dist
,
and niter
. proj
is the projected point, dist
the
distance from the subspace (and niter
the number of iterations
if positivity of the coordinates was requested.).
Some timings show that these implementations are to a certain extent competitive with direct applications of quadprog.
Hans W. Borchers, partly based on code snippets by Ravi Varadhan.
G. Strang (2006). Linear Algebra and Its Applications. Fourth Edition, Cengage Learning, Boston, MA.
nullspace
, orth
#-- Linear projection --------------------------------------------------
# Projection onto the line (1,1,1) in R^3
A <- matrix(c(1,1,1), 3, 1)
B <- matrix(c(1,0,0, 1,2,3, -1,0,1), 3, 3)
S <- linearproj(A, B)
## S$Q
## [,1] [,2] [,3]
## [1,] 0.3333333 2 0
## [2,] 0.3333333 2 0
## [3,] 0.3333333 2 0
# Fourier expansion': sum(<x0, a_i> a_i /<a_i, a_i>), a_i = A[ ,i]
dot(c(1,2,3), A) * A / dot(A, A) # A has only one column
#-- Affine projection --------------------------------------------------
# Projection onto the (hyper-)surface x+y+z = 1 in R^3
A <- t(A); b <- 1
x0 <- c(1,2,3)
affineproj(x0, A, b) # (-2/3, 1/3, 4/3)
# Linear translation: Let S be the linear subspace and A the parallel
# affine subspace of A x = b, a the solution of the linear system, then
# proj_A(x) = a + proj_S(x-a)
a <- qr.solve(A, b)
A0 <- nullspace(A)
xp <- c(a + linearproj(A0, x0 - a)$Q)
## [1] -0.6666667 0.3333333 1.3333333
#-- Projection with positivity ----------------------- 24 ms -- 1.3 s --
s <- affineproj(x0, A, b, unbound = FALSE)
zapsmall(s$proj) # [1] 0 0 1
## $x : 0.000000e+00 3.833092e-17 1.000000e+00
## $niter : 35
#-- Extended Example ------------------------------------------ 80 ms --
## Not run:
set.seed(65537)
n = 1000; m = 100 # dimension, codimension
x0 <- rep(0, n) # project (0, ..., 0)
A <- matrix(runif(m*n), nrow = m) # 100 x 1000
b <- rep(1, m) # A x = b, linear system
a <- qr.solve(A, b) # A a = b, LS solution
A0 <- nullspace(A) # 1000 x 900, base of <A>
xp <- a+drop(A0 %*% dot(x0-a, A0)) # projection
Norm(xp - x0) # [1] 0.06597077
## End(Not run)
#-- Solution with quadprog ------------------------------------ 40 ms --
# D <- diag(1, n) # quadratic form
# A1 <- rbind(A, diag(1, n)) # A x = b and
# b1 <- c(b, rep(0, n)) # x >= 0
# n <- nrow(A)
# sol = quadprog::solve.QP(D, x0, t(A1), b1, meq = n)
# xp <- sol$solution
#-- Solution with CVXR ---------------------------------------- 50 ms --
# library(CVXR)
# x = Variable(n) # n decision variables
# objective = Minimize(p_norm(x0 - x)) # min! || p0 - x ||
# constraint = list(A %*% x == b, x >= 0) # A x = b, x >= 0
# problem = Problem(objective, constraint)
# solution = solve(problem) # Solver: ECOS
# solution$value #
# xp <- solution$getValue(x) #
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.